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NONLINEAR BEHAVIOR OF ACOUSTIC WAVES 
IN COMBUSTION CHAMBERS 
F. E. C. Culick 

ABSTRA CT 

This report is concerned with the general problem of the nonlinear 
growth and limiting amplitude of acoustic waves in a combustion chamber. 
The analysis is intended to provide a formal framework within which prac- 
tical problems can be treated with a minimum of effort and expense. There 
are broadly three parts. First, the general conservation equations are 
expanded in two small parameters, one characterizing the mean flow field 
and one measuring the amplitud€; of oscillations, and then combined to 
yield a nonlinear inhomogeneous wave equation. Second, the unsteady 
pressure and velocity fields are expressed as syntheses of the normal 
modes of the chamber, but with tmknown time -varying amplitudes. This 
procedure yields a representation of a general unsteady field as a system 
of coupled nonlinear oscillators. Finally, the system of nonlinear equations 
is treated by the method of averaging to produce a set of coupled nonlinear 
first order differential equations for the amplitudes and phases of the 
modes. These must be solved numerically, but results can be obtained 
quite inexpensively. 

Subject to the approximations used, the analysis is applicable to 
any combustion chamber. The most interesting applications are probably 
to solid rockets, liquid rockets, or thrust augmentors on jet engines. 

The discussion of tuis report is oriented owards solid propellant rockets. 
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I. INTRODUCTION 

The purpose of this analysis is to develop a suitable framework for 
studying the growth and limiting amplitude of acoustic waves driven mainly 
by interactions with combustion processes. Although the emphasis here is 
on the problem as it arises in solid propellant rocket motors, other cases 
can be treated in the same way. For example, since sources of mass and 
energy within the volvime are accounted for, unstable waves in liquid rocket 
motors and engine thrust augmentors may be regarded as special cases. 

The principal distinguishing features of the solid propellant motor are the 
source of mass and energy at the boundary, and the non-uniform flow field* 

A primary motivation is to produce analytical results which may be 
used to interpret data. For many practical situations, elaborate numerical 
computations based on the governing differential equations are inappropriate 
owing to uncertainties in the required input information. The essential idea 
pursued here is to convert the governing partial differential equations to a 
set of ordinary nonlinear differential equations in time, for the amplitudes 
of the normal modes of the chamber. The way in which this is done is very 
strongly conditioned by previous work on the linear stability of the normal 
modes [Culick (1973 -1975)] and constitutes a development and extension of 
recent work on nonlinear behavior, Culick (1971). The precursor of this work 
was based on the observation that an oscillatory motion in a solid propellant 
motor very often exhibits quite clean sinusoidal behavior even when the 
amplitude attains a limiting value much larger than those at which nonlinear 
effects are clearly evident in, for example, acoustic resonance tubes driven 
at room temperature. This suggested that the acoustic field might be repre- 
sented approximately in the form of a standing wave having time -dependent 
amplitude. 
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^ « T)(t)vKO (1.1) 

u' « V4i(?) (1.2) 

Yk^ 

where k = w/jT is the wavenumber. The true wave structure in space is 
distorted by fractional amounts of the order of the Mach number of the mean 
flow, which, as shown by Culick (1971), need not be explicitly determined 
within the approximations used. The problem comes down to finding the 
amplitude Ti(t), Thus, the unstable wave is regarded as one having a fixed 
shape in space, but the amplitude varies in time. For example, if one 
exairines the fxmdamental mode of a T-burner, i|< = cos (z/L) and the mid- 
plane of the chamber is always a nodal plane. 

The argument in the earlier work led to the nonlinear equation for n, 

T| + + iif(Ti, -n) = 0 (1.3) 

where, partly by assumption and partly from the analysis 

f(Ti,i'i) - -2a + pJtiI + +Yjliil + Y2lill^ (1.4) 

Equation (1.3) describes, as one would anticipate on physical grovuids, a 
nonlinear oscillator, and it is often possible to attach a physical interpre- 
tation to the coefficients a, p^, p^. Yj, Y 2 * 

Approximate solutions to Equation (1.3) have been constructed both 
by the method of averaging, Krylov and Bogoliubov (1947), Bogoliubov and 
Mitropolsky (1961 ), and by expansion in two time variables, Kevorkian (1966). 
Although differing in certain details - for example, higher approximations 
are more easily constructed by using two time variables - the results 
obtained by the two methods are equivalent. In either case, the first 
approximation has the form 


Ti(t) = £7(t) sin(wt + cp (t)) 


(1.5) 



-4- 


The amplitude (7{t) exhibits the correct gross behavior; during growth, (7{t) 
increases from an arbitrarily small value, progressing through a period of 
linear behavior (7 ^ exp (at)^, and ultimately leveling off at some limiting 
value determined by the coefficients in f(r], r \ ). 

However, the approach just described fails in what appears to be an 
important respect. That is, even to second order in some small parameter 
characterizing f, no even harmonics are generated. This is not only con- 
trary to observations in solid propellant rockets, and expecially in T-burners, 
but it cannot be correct if nonlinear effects associated with convection (i. e., 
those represented mainly by the term u • ^u ) are present. 

It is therefore necessary to co’^struct a new analysis. Because the 
results based on the simple analysis just described do in fact exhibit some 
important features of the behavior, it is reasonable to examine modifications 
and extensions of that approach. 

The basic idea here is to permit explicitly, from the beginning, the 
presence of all possible standing waves. This really amounts to stating 
that an arbitrary unsteady field can be synthesized of its Fourier components. 

xt> 

Equations (1.1) and (1. 2) are replaced by the expansions 

oo 

= E ri.(t)4i.(r) (1.6) 

i=0 


S' = E 


-Hid) 


i=l Tk. 


Vilj^(r) 


(1.7) 


The total pressure density and velocity fields are of course average plus 
fluctuation fields 


p 

= p + ep' 

(1.8) 

p 

= p + e p' 

(1.9) 

u 

= + eu' 

(1.10) 


The term i=0 in (L 6) is simply v. presenting a shift of the average pres- 

sure. There is no corresponding velocity fluctuation, so a term i=0 does not 
appear in (1. Only in 510 will the influence of Hq ^ 0 accounted for. 
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where fx is introduced as a dimensionless quantity measuring the magni- 
tude of the near flow speed, and e is similarly a measure of the amplitude 
of the oscillations. Both fx and e are essentially parameters for bookkeep- 
ing (see § 2). 


In § 2 the procedure for constructing the nonlinear wave equation 
is outlined. Mainly the three-dimensional problem is discussed, but some 
contributions arising from the corresponding one -dimensional analysis will 
be incorporated. There are two features distinguishing this analysis from 
previous works treating nonlinear motions in liquid propellant rockets: 

sources of mass, momentum, and energy at the burning surfaces are inclu- 

* 

ded; and the mean flow field is non-uniform. 

The expansions (1 . 6) and (1 . 7) are introduced in the nonlinear wave 
equation, and in § 3 a set of equations for the time -dependent coefficients 
T]n(f)« Each equation of this system represents the motion of a simple 
forced oscillator, where the force depends both linearly and nonlinearly 
on all the 


d Ti 


dt 


^ 


( 1 . 11 ) 


If only the linear terms are retained, then one can extract from (1.11) all 
known results for linear stability analysis. The nonlinear terms arise from 
the gasdynamics in the chamoer, the combustion, and other processes. 

Only the contribution from the gasdynamics can presently be given simple 
explicit forms. Owing to the manner in which the problem has been formu- 
lated here, many of the nonlinear terms represent coupling between the 
modes. It should be noted also that terms representing linear coupling arise 

as well, both from the gasdynamics and from the combustion processes. 

The average pressure, density, and temperature are assumed to be uniform, a 
realistic approximation to the situation in rocket motors. To treat certain types 
of thrust augmentors, one must account for nonuniform average temperatures 
and densities, which ran be done within the framework developed in this report. 
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In §4 an approximate means of solving the set (1,11) is discussed. 
While it is true that if is given, the coupled equations can be solved 
numerically, this may be a relatively expensive procedure. The purpose 
here is to provide a considerably faster and cheaper means of obtaining 
the informucion desired. Here the technique used is essentially that termed 
generically the "method of averaging. " It is based on the assumption -- 
almost always valid for the unstable motions encountered in practice- - that 
the motions exhibit relatively slowly varying amplitude and phase. Thus, 
the functions represented as 

T) (t) = (7 (t)sin(co t + cp (t)) =A (t)sin 03 t + B (t)coso3 t (1.12) 

'n n ^ n n' ” n' ' n n' n ' ' 

According to the basic assumption used, the quantities ^7 , <¥> , A , B 

n n n n 

suffer only small fractional changes during one period of the oscillation. 

The analysis then produces coupled first order ordinary differential equa- 
tions, a system which is cheaper to solve than the system of second order 
equations. 

The system ol first order oquations is valid for problems in which 
the frequencies of the higher mo^cs are not necessarily integral multiples 
of the fundamental frequency. In §4 the general equations are given. As 
an elementary example, the motions of two coupled pendula is analyzed in §5* 
this shows the familiar beating of the oscillations, a feature which seems 
not to have been accommodated by previous applicavions of the method of 
averaging. 

Many practical problems involve purely longitudinal ('^organ pipe'’) 
modes, for which the frequencies are integral multiples of the fundamental. 
The system of nonlinear first order equations simplifies considerably for 
this case, treated in §6. For applications, it is necessary to incorporate 
representations of processes responsible for the loss and gain of energy 
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by the waves. An approximation to the interactions between pressure waves 
and surface combustion is described in §7. One way of handling the loss of 
energy due to particles suspended in the gas is covered in §8; the results 
show favorable comparison with more exact numerical calculations reported 
elsewhere. In §9, linear and nonlinear viscous losses on an inert surface 
are examired. Associated with the nonlinear unsteady motions there is also 
a change in the average pressure; this is discussed in ^10. 

Several examples of unstable motions in motors and T-burners are 
covered in 511. The cases have been chosen for comparison with numerical 
results previously reported; again^ the agreement appears to be quite good. 

As noted above, the analysis has been strongly motivated by previous 
work on the linear stability of motions. In §12 the connection is discussed. 
One of the attractive features of the formulation cf the noxilinear behavior 
is that the more familiar linear results are not merely accommodated, but 
explicitly incorporated and used. An interesting and important unsolved 
problem concerns the influence of coefficients characterizing linear behavior 
on the nonlinear behavior. The coefficients are proportional to the real and 
imaginary parts of the complex wavenumber computed in the linear analysis. 
Nonlinear behavior appears to be quite sensitive to their values; a few ex- 
amples are included in the brief discussion given in §12. 2. 

Analysis based on expansion in normal modes with time -dependent 
coefficients have earlier been reported for unsteady motions in liquid pro- 
pellant rocket motors [e. g. Zinn and Powell (1970), Lores and Zinn (1973) 
and other works cited there]. Results were obtained for specific problems 
by solving the second order equations for the amplitudes. Reduction to a 
set of first order equations was not effected. Thus, the computational costs 
must be substantially greater. Moreover, interpretation of the formal rep- 
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sentation, and incorporation of processes such as particle damping and sur- 
face heat losses appears to be somewhat more difficult than for the analysis 
developed here. It is likely that the techniques discussed in the references 
cited above and those discussed in this report should produce the same sec- 
ond order equations for the same problem. This has not been verified, 
n. CONSTRUCTION OF THE NONLINEAR WAVE EQUATION . 

The nonlinear wave equation for the pressure is constructed by 
suitably combining the conservation equations and the equation of state. 

For applications to solid propellant rocket motors, it is necessary to 
treat the mediiim in iiie chamber as a two-phase mixture of gas and parti- 
cles. Culick (1974) has outlined the steps necessary to produce the equa- 
tions for the velocity and pressure 

^4 

p If + pu . Vu + Vp = 6 Fp - a (2. 1 ) 



u • Vp + Y pV - u 


k • 


F + u 


+ (e - e )w 
' po o' p 


+ (Q + 6Q ) + ( 1 + K)C Tw 3 
P ^ P 


( 2 . 2 ) 


Here, p is the density of the mixture, p = p + p p and h is the ratio of the 

P 8 

mass of particulate i latter to the mass of gas in a unit voltime of chamber: 

H =: p / p . It will be assumed throughout that K is a constant in both space 
P 8 

and time. With C the specific heat of the particulate material, the proper- 
ties of the mixture are: 


C + KC 
C = -I 

V 1 + K 


p = Pg(i + ^ > 


c = 

p 


C + KC 
-E 

1 + H 


c 

7 = -=^ 

c.. 


(2.3) 


= yR T = ^ 
P 




-9- 


The equation of state is 


p = RpT = 


Rp T 
1 + K 


(2.4) 


where R is the gas constant for the gas only, and R = = R/(l 

Note that ( ) is used In (2. 3) and (2.4) to denote certain material proper- 
ties of the mixture. Later the same notation will be used to denote time 
averages of variables. 

The force of interaction and heat transfer between the gas and 
particles is represented by 


F = 
P 




3e 

, _JE 

p 8t 


9u ^ 

P “37^ + p u • 
^p ct ^p p 

%] 


(2. 5) 

p u • « 

-p c 

^p 

r 3 t ^ 1 

(2.6) 


The differences between the local values of velocity and temperature are 

T (2.7) 


6a = u - u 
P P 


6T = T 
P P 


Then the differential force and heat transfer are 

*^p = 'Pp[ ^ <"p * % - (2- 8) 

r36T 1 

«Qp = -•'pC[-5r* + <v " • '■^’J ’> 

All symbols are defined in the list of the end of the report; additional details 
leading to the forms quoted here may be found in Culick (1974). 

The nonlinear wave equation is found from (2. 1 ) and (2. 2) as 

.2_ 


3t" 


- Yp7 • ^ - Yp7 • (u • Vu) u+-^(u* 7p)J - Vp7*^-^-^ 


3P 


( 2 . 10 ) 


where 

The corresponding results for purely one -dimensional problems Lsee Culick 
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(1973, 1974) for further details of the formulation] is 




c 




- YP 


S 3z\®c p ) 


8P, as, 

* -dT ~sr 


with 


Xj = (U-U )w + |-[ujm^^dq + u J m^^P^dq] 


( 2 . 12 ) 


(2.13) 


P, =^-r(u - u)F +uX, + (e -e )w +(Q+6Q ) 

1 C L P P 1 ' po o' p ' ^p' 


+ (1 + K)C^TWp] 


(2. 14) 


h = k)yR(T+ AT) + i ^ 

^ r|=^T +4r(..pb" - »p'] Pp^V’ 


(2s 15) 


It has been assumed, to obtain the form shown for Sj»that AT, the difference 
between the temperature of the material at the edge of the combustion zone 
and the average value in the chamber, is the same for both gas and particles. 
This is not an essential assumption, but is done here to simplify the formulas 
somewhat. Note that except for extra terms in Xj and Sy there is a one- 
to-one correspondence between the terms of (2. 10) and (2. 12). 


There are two ways of proceeding towards soluble problems. A 
formal expansion procedure can be applied directly to the complete wave 
equations, (2.10) and (2.11); or the first order equations (2.1) and (2.2) 
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can be expanded, and the wave equation formed. The results are identical, 
but because the second calculation is somewhat simpler, it is given here. 

Si-bstitute the expansions (1.8) - (1.10) into (2. 1 ) and retain only 
terms to second order in e to find: 


It + = <^p' = -l^Lu.Vu'+u'.Vu] - e[u'-Vu'+ ^ 3 +7^(6F'-a') (2.16) 

dr p p dt ®^p p 

-|B + YpV* u' = “p[u • 7p* V* u ] - €[Yp>7 *u* +u** 7p*] + P* (2, 17) 


where P* is the fluctuation of P. To this point terms to all orders of p have 

hficn retc-.ined, and no assumption has been made about the ordering of 

o' and P*. Eventually only terms to first order in the mean flow speed will 

be retained^ so the mean pressure and density will be constant. For 

❖ _ — «2 / 

simplicity, use that fact now and write p = p^, p = p^, a = ^Pq/Pq- 

The nonlinear wave equation may be formed now by differentiating 
(2. 17) with respect to time, and substituting (2. 16) into the second term 
on the eft hand side. Some rearrangement gives 




-fx|po^ • (u • 7u' + u'« Vu ) U* 7^ - ^ 7^ 7* ul 
' ® a a ' 

- {po^‘ («'• ^ It <p'^* «' ) - Z2 ^ • ^p' > + <p'lr' )} 




(2.18) 


The boxindary coIK^’~ion accompanying (2. 18) is found by taking the component 
of (2. 16) norr^.al to tlie boundary: 


n • ^p' p^(u * 7u' +u'* 7u)* n - e(p^u'* 7u' + p' ^-)* n 


(2.19) 


+ -7 (6F ' -a‘)*n 
e ' p ' 

The cne -dimensional counterparts of (2. 17) and (2. 18) are easily 
deduceu by replacing 7 by 9/az, and the divergence 7*v of a vector v by 
Note that p stands or the average value of the flensity of the mixture : 

P =p +p -p (Hk). 

o ' po ^ go ' go 
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S ~^8/8z(S V ). Also, o is replaced by S, and P by P, +S, . 

For computations of linear stability .only terms of order p appear on 
die right hand sides cf both (2. 18) and (2. 19). The solution for the n^ 
mode then has the form (1.1), with r|(t) a simple exponential, 

at _ 

p^' = p^ e “ t|ii^(r) (2.20) 

Thus, the first influence of the perturbation produces a growth or decay; 
a^ is proportional to the average flow speed, characterized by p. There is 
a small change of the frequency, but the spatial structure of the mode, 
is undistorted in first approximation. Calculations of-that-sort can be 
extended: the next approximation provides a formula for a^ correct to 

second order in the mean flow speed, and the first order distortion of the 
mode shape. 

The purpose here is to develop a description of nonlinear temporal 
behavior associated with second order acoustics. The perturbations associ- 
ated with combustion and mean flow are the same as those accounted for in 
the treatment of linear stability. In order that only linear terms in p should 
appear in (2. 18) and (2. 19). it is necessary that the undistorted mode 

shapes should be used on the right hand sides. Formally, this step implies 

2 

that terms of order p , ep and higher are neglected compared with those 
of order p and e , TMs corresponds to the following limit process applied 
to the small parameter p and e. 

Consider, for example, only the two terms on the right hand side of 

(2.18), 

p I p^V • (u • 7 u' )| + e jp^7 • (u‘ • 7u' )| 

The acoustic velocity field with first order distortion is 

a' = u ' + pu, ' 
a '^1 
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where u ' is the classical acoustic field. Substitution gives 

St 

[IpqV . (5 . Vu^' ) J . 7u^‘ )| 

p^V • (u • 7Uj ' )| + e |p^7 ♦ (u^' • 7uj ' + Uj ' • 

+ Jpo7-(Si' •7ui')jj 

Now let p, € -• 0 but e/p.~ 0 (1). Only the first terms survive, containing 
the undistorted parts of the acoustic field. Similar reasoning applies to 
the other terms on the right hand sides of (2. 18) and (2. 19). The procedure 
can be extended to hi-^her order in both p and e. 

The expansions (1.6) and (1.7) are therefore appropriate. How 
they are used to give a coupled set of equations for the amplitudes T]^(t) 
is described in the next section. 


m. ORDINARY NONLINEAR EQUATIONS FOR THE AMPLITUDES 
3. 1 Construction of the Equations 
Define the functions 

h = - p 7* (u • 7u' + u'* 7u) + ■^u* 7^- + — - 7*u 

p o —2 ot —2 ot 

a a 

*(u' -VuM + S. 4(p'V-u')+-i= -|-(u'-7p')-7-(o'l^'l 


h = -p V* (u' 




, 8u' . 

= Po'5t 


f = p (u • Vu' + u'* 7u) • n 

f = . i (6F ' . 0')- n 
v e p 


(3.1) 


e -o ' ^ at at ^ 


(3.3) 


(3.4) 

(3.5) 

(3.6) 

(3.7) 

Then the nonlinear wave equation (2. 18) and boundary condition (2. 19) become* 


The ordering parameters are hereafter suppressed except in h^ and f^. 
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a dt 


(h +h +h ) 
U e V 


(3.8) 


n . 7p' = -f - f - - f 

s p e V 


(3.9) 


Multiply (3. 8) by the mode shape for the nth mode and integrate over 
the volume of the chamber. The first term on the left hand side can be 
re-written using Green's theorem, and (3. 9) is tlien substituted. Some use 
must be made of the following properties of the 


+ h = 0 
^n n ^n 


(3.10) 


n • = 0 


rd;.4, dV = 6 . E ^ 
" u^n ni n 


If 2 2/^ 

K = o> / a 
n n ' 


These operations lead to 

»2.,i 


at 


a^0ili [f +f + f -+ f 3 dS 

^ ^n s p s V 


(3.11) 

(3.12) 

(3.13) 


(3.14) 


Because of tlie orthogonality (3. 12) of the substitution of (1 . 7) 
in the left hand side jf (3. 14) gives 

2r . 2 T —2 


p E ‘^[tL+oj ti 3 = - a^ filj [h +h-+h 3dV 
n n n n J^n'’ u. e y 


(3.15) 


- 2 , 
- a 


. [f +f + f +f ]dS 
n s |i. e V 


where 


For the one -dimensional formulation, the equation corresponding to (3, 15) 
is 
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(3.17) 


where 


lli 


‘le 


= J 

o 

(3.18) 

a 8 u 8^p' .7 8p* d . 

a c 

(3.19) 

,8u\ , Y 3/.9/P. IV. 19, ,9p\ 1 8 

8t <P 8z<V ) 2 .-^ 

a S a a S 

c c 

(ScP't-') 

(3.20) 

•?{r2 

a c 

(3.21) 

f - ^ 

Is ~ 9t 

(3. 22) 

*lp= f“oK<""'> 

(3. 23) 

^le = Po" ^ + P ^ 

(3.24) 

f, = - i (6F ' - El ) 
Iw s P 1 

(3.25) 


3. 2 Evaluating the Linear Terms 

After some re-arrangement, one can establish the following identities: 

u' )4,^dV - J[S' XVXS]- 7^,_^dV 

+ V^7-3]dV 
a 


(3.26) 
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*^Ort ao c 

(3.27) 


/WV + i 2+nlf ' ■" <‘^p’ ■"' ' • ’•I'n^ 


(3.28) 


— "^1 1 ft 

K’'lw^'' + [ Wlw\ = - T ( +t It l^i +sp +(6F; - g/]S^dz 


o a 


(3.29) 


The terms in (3. 28) and (3. 29) represent the influences of both surface and 
residual combustion within the volume; the contribiitions associated with 
condensed material sut^pended in the gas; and, for the one -dimensional 
problem, the effects o' flow entering fiom the lateral boundary. Later(§3, 3) 
the peculiarly one -dimensional terms will be combined in the three-dimensional 
problem according to the prescription discussed by Culick (1975), To sim- 
plify writing, define the functions H, Hj and Lj, which appear in (3. 34) and (3. 35); 


1 r 1 ^^1 1 

^1- 


9S, 


= kLh ~sr ■ (^1 ■ 


dz 


(3. 30) 

(3.31) 

(3.32) 


The acoustic quantities appearing in (3. 26) and (3. 27) are approxi- 

ji- 

mated by the expansions (1.6) and (1. 7) to give 


p oo Ic^ 

/\l/ hdv + < 5 \)/ f dS =— S i f — ^ t u • V -\|f. u • V \|^ 

•' a Ji u — . ,1 L-l L, 2 ’^n 

^ y 1=1 k. 


dv 


( X 7 xu)\|/^* dv 


+ (y-1) /'t'-t 7*U dv 


1 K 

-2 F 


1 n 


9p' ^ A 

t u . n dS 

dt n 


(3. 33) 


’■'as noted in connection with (1.6), r\^ will be assumed to be zero throughout 
except in 5 10. 
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00 


I "iM m { I [-|^ n “ -5T- - +m“ ■S]®c'*" 

^ m 

o c 


i: 


(3.33) 


The terms containing £ and f, are conveniently combined with the last 

S X s 

terms of (3. 32) and (3. 33). Then with the preceding results. Equations 


(3. 15) and (3. 17) become 

oo 


^n ^*n “n = "I ^idC’^T '^n 'I'i u- V Jdv 

i^l i 

- J(V\|^i X VXu)>|f^ dv 
+ (y-l)Jtj'kj^7-u dv]" 


, d 

P '"bl 

♦nar 

Po“b* p-] 


(3. 34) 


_ JL 


{ I ♦„ '’e + 1 1,, ds} 


._2L 


Jh 


dv 




m^l 


m 




(Continueci) 


The signs of L J the surface integral have been changed so 
uhn = and u*n = U|^ are positive inward* 
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(Continued from Page 16) 


a "‘o 

_ ^ L ^ 

o 

+ r H|S dz “f r Li* S dz 
pJ Ic o ^ I c 


(3. 34) 


3. 3 Surface Terms in the One -Dimensional Problem; Combinir .: 
the One- and Three-Dimensional Problems 
It is the function Lj which contains the terms found in the one -dimensional 
approximation. Only the linear forms will be treated in this work. With the 
definitions (2* 13) and (2. 15), Lj to first order is 

o o 

- [u- J dq + u* J i^^P^dq]-^^ (3. 35) 


It is readily established. Equation (3. 11) of Culick (1974) that 

+ • < 
o o pa 

Then because the integral over dz dq is the integral over the lateral sur- 
face, we can write 


(3. 36) 


dz= xjj tp.^ldS 


|_fj- ^ ^ dS dS 


,_ii_ ^(pids 

p Jtl p dz b 
^ 

* . r T 

The sign of u^ in I J has been changed so that Uj/ is positive inwar<l 
at the end surfaces. 
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[t is a good approximation for most practical probl<-ims that the temperature 
fluctuations within the chamber are nearly isentropic. With this assumption, 
and the approximation that the composition of the material leaving the com- 
bustion zone is the sante as that within the chamber (so p = (1+K)p ), the 

o g 

last identity can be written 


-2L 


J Ljdz = +p':-^)ds 


O O 


/— , t 9 


^o 

dt. _ 


1 ^. 37 ) 




The terms multiplied by (y-1) in (3* 34) and (3. 37) combine to give 

00 L 

w 


O '-*'! VJ 


(3. 38) 


m- 1 


g 


where use has been made of the continuity equation for one -dimensional flow. 


Substitution of (3. 37) and (3, 38) into (3, 34) then leads to 

00 , Z 


‘ -1 ’’ miJtrrh " - ♦ m " 


m=i 


m 


p dt 


° “s 

+[Sch<'>o“b^P':^] } 


- L 


(Continued) 

This assumption which, at the expense of substantially more labor, < an 
be relaxed, means that temperature or entropy wa-^es (convected by the 
average flow) are not represented in the volume of the chamber* Nonis “in- 
tropic temperature fluctuations at the surface are accounted for (AT ^ 0)* 
There is at present no evidence or calculation showing how important this 
inconsistency might be. Essentially v/hat is included is the direct influence 
of nonisentropic surface combustion on the acoustic field, but interactions 
between the entropy Wtives and the acoustic waves are ignored* 



(Continued from p. 18) 
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I ^ 

m=l m 


+ -3i r H, S dz 
p 1 c 

o 


(i. 39) 


Note that the surface integrals in (3. 39) extend along the lateral 
boundary only, and that the velocities u^, appearing in those integrc^ls 
are positive inward. '!^he second set of brackets containing surface terms 
obviously correspond exactly to the surface integral in (3. 34), According 
to the argument proposed by Culick (1974) these should be written in the 
general case as* 

where U|| ' and \ij^ both stand for the formula (3. 36). 

Moreover, the surface terms contained in the last brackets of (3. 39) ccn be 
incorporated in the three-dimensional problems in the form 
^ T] 

I ri S ♦n* V® + ^ ^ "b - • 


Consequently, with pr< per interpretation, all problems are represente ! by 
the equation 



n 







V\l 

n 



dS 


+ ^jHdv + -^H Si 

^O 


(3.40) 


The symbols u*, u*|j, 6 , 6 y are defined and discussed in § 4 of Culick 
(197 5). 



where 
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D . 
n ni 


rr^n ^ 1 

1 -y i u . V 4 - - \|f u • V i dv 
l_ ^ ^ n 1 1 n j 

i _ 

+ (y - 1 ) / + "T ^ ^ 1'*^ 


(3.41) 


R = (u'j^6x+ a'||6|i) 


(3.42) 


That the nonlinear .erms for the one- and three-dimensional problems corre- 
spond exactly, and therefore can be accommodated as sho^im in (3.40), is 
demonstrated in the next section. 

3.4 Nonlinear Terms 

To the order considered here, there are no nonlinear terms explicitly 
dependent on the mean flow. The representations will therefore be directly 
useful for the classical problem of nonlinear waves in a resonant tube. No 
special considerations are required for the one -dimensional problem. With 
the definitions (3.2) and (3.6), 


I h„dv + ff\jr fdsj 

e D i ^ m £ Js ^n £ J 


= _r 


• P —* -♦ 

JV it * 0 u*eVu’ + p* -57 — dv 
n L o ^ at _ 


(3.43) 




It is within the approximations already introduced to use the zero order 
acoustic approximations in these integrals: 


^ - 
at 


-VP V* u* 


(3.44) 
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Then (3. 43) may be written as 

Ig = y (u'*Vu') - ^j^{y(V-u')^+u'.V(V.u')}]dv 

- -r^ f (p' Vp') + +^C(Vp')^ +7 p' 7^P'}] dv (3. 45) 


—2 

In the first term of the second integral^ the approximation p* a p* has 

♦ 

been used; this is consistent with (3# 44). Now substitute the expansions 
(1.6) and (1. 7) to find, with the terms involving r\^ dropped. 


OC 00 


) y Fa .. •n.n. +B ..-n.Ti.l 

e n L. Li L nij r j nij i jj 


i=l i=l 


(3.46) 


where 


“J yE^ k.^k.^ " " J 

' -n t 1 


n 1 j 


2 2 2 

- I' [y k. k. 'if.i/. - k. ]dv 


-2 


B .. 
nij 


-k, i.t.}]dv 


= ^ /[7^ + ')„{(vxk.)*(^'l'.)- y j Vj 

yE^ n 1 J n 1 J 


(3.47) 


(3. 48) 


Four integrals appear in these definitions: 


II 

/ 7^^* (V\|^.-7(7\|f^)]dv 

(3.49) 

II 

J * (\|(. 7\|r. )dv 

'^1 j 

(3.50) 

-ij _ 

"3n 

f yh 7 >!(.• 7>i(. dv 
1 J 

(3.51) 

I .. = 
niJ 

/ >1^ >1^. t- 

n 1 j 

(3.52) 


More detailed consideration has been given by Chester (1961) to the use 
of acoustic approximations in the representation of finite amplitude waA es. 
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Then (3. 47) and (3, 48) are 


— 2 2 2 
y E k. k. A . - 
' n 1 j mj 

ye^ 

n B .. 

-2 mj 


= + k^ I - y k? k^' I 

In j 3n ^ 1 j r 


+ I, - y k. I .. 
2n 3 j mj 


(3.53) 


(3.54) 


One can eventually reduce (3.49)-(3. 51) to multiples of 

xij 1 ri.4 n 2 . 2 . 2 ,^ 1 , 2 „ 2^,2 , 2 .. ^ 1 „ 2 , 2 ... 2 ^, 2 , 2 ,_ 

I,'' = xLk. -(k -k. ) JI .. = (k. +k. -k )I .. + x(^- “k- )(k- )I •• 

In 4 J n 1 ' mj 4 n ' i j n mj 4'j i''i j n' mj 


jj 1 „ 2 , 2^, 2„ 1 1 2 _ J. 1/. 2 , 2„ 

IJ = Tj- (k. -k. +k )I . . = X k I . . + •^k. -k. )I . . 

2n 2 j X n n’j 2 n mj 2 j i mj 


1 , 2 


1/, 2 , 2, 


(3.55) 

(3.56) 


= x(k-^ + - k^)I .. 

3n 2 1 j n mj 


(3.57) 


The secon . parts of (3- 55) and (3. 56) show and decomposed intc pieces 
which are respectively symmetric and antisymmetric under interchange of 
the indices (i, j). With these results, (3.53) and (3.54) become 


A .. = ^ [(k.^+k.^)^ - k^ - 4 y k.“k.^] 

1 J n 

+ -y — J (k -k.^ ) (k ^ +k.^ -k ^ ^ 

2yk.^K^ J ^ J ^ " 

1 j n 

B .. = 5?3J. (k. +k ) t ySH (kf - fc. ) 

2yE^^ J ' 2yE_,^ J ‘ 

Again, the two parts on the right hand sides are respectively 'ymnetric 


(3.58) 


(3.59) 


and antisymmetric in (i, j). In ^ 4, Equations (4. 32) and (4. 33), the following 


combinations will arise: 
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A .. to), a). ± B .. = ^ I 

nij X J mj 2yE^^ 

-2t 

a I 


a^I . . (k.‘'+k.“)“-k 

I 1 2 


2., 2,2 , 4 
n 


2 k. k. 
X J 


- 2yk.k.±(y-l)(k^+k.^)] 


,-k.^+k.^-k^ 


€3. X - o r n 


^ J 


(3.60) 


It should be noted that there is yet no restriction to a particular 

ceometry; the form of the chamber influences the numerical results through 

the values of the k. and the integral I ... Because some terms in I , eq. 

1 ® mj e 

(3.45), contain p', there are non-zero values for i = 0 or j = 0. These are 

associated with the DC shift of mean pressure due to nonlinear processes; 

this is treated separately in § 10. 

Finally, because the nonlinear terms (3.46) are summed over all 

(i» j). 2 tnd the coefficiei ts ^nij multiplied by functions which are 

symmetric in (i, j), only the terms containing the symmetric parts of 

3 . will survive. Thus, only the first brackets in (3. 60) will be required 

niJ 

Jor later calculations. 


IV. APPLICATION OF THE METHOD OF AVERAGING 

Most of the ter.ms on the right hand side are such that (3. 40) may be 
brought to the form 


^ 2 
T) + CO 
n n 


n 


F 

n 


(4.1) 


vith 


00 


00 00 


n 


- ) CD . Tj. + E , Tl. ] - ) y [a . . . .T).T1. ] 

L ^ m 1 ni 1 L L njj 1 j nij i j 

i=l i=lj=l 


(4.2) 


Other contributions (for example proportional to It). I, , etc. ) may 

arise; they are easily handled within the framework to be d.jscribed now and 
need not be considered explicitly. 



-25- 


The set of coupled equations (4. 1) can be solved numerically, but at 

considerable expense. It is the intent here to reduce the second order 

equations to first order equations, for which solutions may be calculated 

quite cheaply. The basis for the approach taken is the fact that 

many of the observed instabilities are essentially periodic with amplitudes 

slowly changing in time. Each mode may thti *fore be reasonably represented 

by the form (1. 12) with Q (t), cp (t), A (t) and B (t) slowly varying junctions 
* nnn xi ^ * 

of time, equations for the amplitudes and phases are found by averaging 
{4, 1) over an interval T which will be defined later. 

Accounts of the method of averaging have been given by Krylov and 
Bogoliubov (1947) and Bogoliubov and Mitropolsky (1961). The results are 
restricted by the condition that must be periodic. This is true for the 
problems treated here if the modal frequencies are integral multiples of the 
fundamental, a condition which is satisfied only by purely longitudinal modes. 
Moreover, solutions are required here for a time interval longer than that 
for which the more familiar results are valid. Both difficulties are overcome- - 
to some approximation not clarified or examined here --by the following 
heuristic development. 

Equation (4. 1) represents the behavior of a forced oscillator, for 
which the motion is 

T) (t) ^ d (t)sin(co t + Cp (t)') = A sino) t + B cos O) t (4. 3) 

n n ^nn^n n n n 

and 

a = [A^ r (4.4a) 

nnn ' 

cp = arctan [B /A ] 
n n n 


(4. 4b) 
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The energy of the oscillator is 



1 2 2 
2 



(4.5) 


Because the oscillator has instantaneous velocity the rate at which work 
is done on the oscillator is ^ F^. The values time -averaged over flie interval 
T at time t are 

t+T t+T 

<«?„> = «V=7K*'n*’ 


Conservation of energy for the averaged motion requires that the rate of 
change of time -averaged energy of the oscillator equal the time-averaged 
rate of work done; 




<^^n> 


(4.7) 


In all that follows the essential assumption is used that the fractional 
changes of the amplitude and phase are small during the interval of averaging# 
The changes in time t are approximately Sl^ and so the asstimption 

is 

T << X , cp T < < 2ir (4. 8) 

n n 

According to (4. 2), the velocity of the oscillator is 


Ti = ui a cos(oj t-np )+ [(? (3 cos(tot+cp )+i? sin(o) t + \|7 )] 
n nn nn nn nnn nn 

The inequalities (4, 8) imply that the terms in brackets are negligible com- 
pared with the first term; the stronger condition is set [see Krylov and 
Bogoliubov (1947), p, 10 ] that the combination vanishes exactly: 


cp cos (OJ t + cp ) + ^ sin(to t ^ cp ) = 0 
n n n n n n n 


(4. 9) 
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Thus the velocity and energy are 

f\ = u (7 cos((0 t + cp ) 
n n n n 


(4. 10) 





(4.11) 


A second consequence of (4. 8) is that when the integrals in (4.6) are done, 

(7^ and are taken to be constant --i. e. they do not vary significantly during 
the interval of averaging. Equation (4,7) therefore becomes 


da 


t+T 


UUf « M 

- -4^ r 


-nr- - ' F cos(a) t'+cp )dt' 

dt W T 0 n ' n n' 


(4.12) 


n 


An equation relating and (5"^ is found by substituting (4, 2), (4. 9), 


and (4. 10) into (4. 1). The time-averaged result is. 


dcp 

dF 


t+T 


CO TC7 
n n 



sin(to t* + cp )dt* 
n n 


(4.13) 


Although a^, cp^ are approximately constant over one cycle, they may vary 
substantially over long periods of time. Equations (4, 12) and (4. 13) are 
then awkward to use. The difficulty is avoided by solving (4, 3), (4,4), 

(4. 12) and (4, 13) for and It is this pair of equations which will be 

used as the basis for subsequent work: 


dA 

n 

dt 

dB 

n 

dt 


t + T 

= — I F cos (0 t' dt' 
(0 T J n n 

” t 

t + T 

= r F sin (0 t' dt' 

(0 T J n n 


(4. 14) 

(4.15) 


For F^ given by (4.2), the first order equations have the form 
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dA / qA \ / dA \ 

r = (irL 


“ar 

dB 


linear 


nonlinear 


«r>_ / dB \ / ax> v 

“ar = (" ar ) +(“ ar ) 

linear nonlinear 


dB 


( 4 . 16 ) 
( 4 . 17 ) 


The linear contributions are: 

E 

nn' n *• "(« 




1 A ^ Ml -Q 


linear 


^ H ^ 'f(W*D .[(f .+h .)A.*(g 
n 2 ( 1 ) ^ 1 ni m i ni 

n 


/dB X 

(“ ar ).. 

linear 


\e .[(g ..h jA.l} 

2 ( 1 ) I ni ^ni ni i ni ni i 
n 

= -|D B +|-iHlA (u).D .C(f ,-h -)B.+(g .~l .)A.l} 

^nnn^u) n2(i) li nx m ni i ni 

n n 

+ 7^^ {e .[(g .-I ,)B.-(f .-h .)A.]} 

2«)_ ni '’ni ni i ni m i 


( 4 . 18 ) 


where 


sin(o}.-fu) )^ 


sin(a).+co )j- 

*»i = ~Trfr- >inC(Mj+u^)(t+j-)] 




sin(o}. ~UJ ) j- 


I . = 
m 


(w.-co ) =- 

' 1 n 2 

sin(o}. -w ) y 
1 n 2 


(COi-W^) J 


sin[((0j-W^)(t+2 )] 


( 4 . 19 ) 


( 4 . 20 ) 


( 4 . 21 ) 


( 4 . 22 ) 


( 4 . 23 ) 


The nonlinear contributions are 

.dA 


( UJX \ 

V ") 


nonline a: 


i=l i=l 

+ D ..lb.. T™^ + d.. T^^l} 

nijLiJ 1“ ij 3- JJ 


( 4 . 24 ) 
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where 


/dB 


00 00 


Q) .. = 5^ I Z {Crfjtij T?.' - Sj ] 

nonlinear n ._j ^ 




■ «1J _ 
1± ■ 


Sin (to +00.) 5* * 

^ cos [{oo +to . )(t + 2 ) 

sin(to -to.)^ _ 

+ K £ cos[<to -to^)(t + 1 )] 

K'"±^2 


(4. 26) 


™nij 

^2± 


sin(oo +W,)^ 

2 ^ sin[ (to^+ to^)(t + 1 )] 


(tOn+tO^)2 
sin(to^-tO^)J 


K‘"±)2 


sin [ (to^ - to ^)(t + j )] 


(4.27) 


T^j 

3± 


sin (to +toJ 5- 

— — sin[(tO^+tO^)(t+^)] 

(tO^+tO^)2 


sin(to -tO,)y 

s ^ sin[(« -w . )(t +i )] 


<“n -“±>2 


1 / 

n ± 


(4. 28) 


T^J 

4± 


sin(to fco )y 

;:p cos[(tO^+tO^)(t + j )] 




sin(c 0 -to.) 5 - 

P - - y - cos[(t 0 ^-t 0 ,)(t + 1 )] 

r 


(4. 29) 


to. = to. + to. 

+ 1 J 


t0_ = (Oj - tOj 


(4. 30) 



b.. 


(4. 31)a, b. c.d 




C-. = T (A.B + A.B. 

ij 2 ' 1 j j 1 


d.. = 4- (A.B. - A.B.) 

ij 2 ' 1 j j i' 


C .. = A .. co.co. - B .. 
mj nij I J mj 


D .. = A .. CO.CO. + B .. 
nij nij 1 J nij 


(4. 32) 
(4. 33) 


The coefficients C .. and D .. are calculated with Equation (3.60); as noted 
mj mj ’ 

at the end of § 3, only the symmetric parts are required. 

These formulas are valid for any geometry; the modal frequencies 

may have any values. Considerable simplification may accompany special 

cases. Most of the following discussion will be concerned with problems 

.nvolving purely longitudinal modes, for which - nWj. 

The interval of averaging remains unspecifi-id. Two possible ties 

are fairly obvious: t = T^, the period of the fxmdamental oscillation; and 

th 

T = T , the period of the n mode. In the first case, each equation is 
n 

averaged over the same interval, while if T = T^, each equation is averaged 

aver its own period. Partly because the argument leading to (4. 14) and 

(4. 15) is not rigorous, there is no wholly satisfyinj^ reason for choosing 

one or the other alternative. If the same interval, T = T^, is used for all 

aquations, then it is necessary, for (4.8) to be satisfied, that the amplitude 
th 

and phase of the n oscillation not change much in roughly n of its own 
periods. This same condition must be met if t = in each of the n equations. 
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For in the equation n = 1, and can be taken outside the integral only if 
they are nearly constant over the interval Tj = n 

In all the problems considered here, the motions are in fact domi- 
nated by the fxmdamental mode: the time scale of the slow changes is 
usually longer than and is essentially the same for all modes. Hence, 
the choice T = is appealing and will be used here. It happens also that 
in some cases the equations are somewhat simpler. Comparison of numer- 
ical results for the two possibilities has not been made. 

V. AN EXAMPLE OF AMPLITUDE MODULATION 
It appears that applications of the method of averaging, and the 
procedure based on expansion in two time variables as well, have been 
restricted to problems for which the function F^, Equation (4. 1 ), is periodic. 
This is true, for example, if the modal frequencies CO^ are integral multiples 
of the fundamental, = nco^ Because the specific examples discussed 
later are, for simplicity, also based on the condition that = n(Oj, it is 
useful to examine a simple case in which the frequencies are not so related. 
This may serve not to prove but to suggest the validity of Equations (4, 14) 
and (4. 15). The practical importance of this conclusion is considerable, 
for one is then in a position to treat nonlinear problems involving tangential 
and mixed tangential /axial modes which are commonly unstable in certain 
kinds of combustion chambers. 

Consider the simple problem of two oscillators linearly coupled and 
described by the equations 

Tl'i +c/jTlj = KTI^ (5.1) 

^2 = *^^1 


(5.2) 
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Equations (4. 18) and (4. 19) reduce to 


dA. 

1 

dt 


^ C(g.. + Jt..)A. + (£..+ h..)Bj 

Zu}^ ij' j ' ij ij' j 


dB. 

1 

dt 


2oJ. 


(f.j - h.j)Aj] 


(5.3) 

(5.4) 


Here, if i = 1, then j = 2 and conversely, giving four equations. For the case 


when the two oscillators have nearly the same frequency, 0 )j-i-a )2 ^ 2(0j, 

CO, -CO- 0, so h. . 1 and f., g., 1.. -* 0. It is then a simple matter to show 

i C IJ 1 1 IJ 

that the A., B^ all satisfy the same equation. 



(5.5) 


The coefficients all oscillate at the "beat frequency," approximately equal 
to K/2c0j. 

For example, if the initial condition is T)j =0, ^2 * then a solution 
is 


Tij = Cj sin( t) sin cOj t 


(5.6) 
(5. 7) 


VI. LONGITUDINAL MODES: CO = nco, 
n 1 

This special case will serve as the basis for many of the examples 
discussed later. The integrands in (4. 14) and (4. 15) are now periodic, with 
period ; the limits on the integrals can therefore be changed from (t, t+j) 

to (0, t). In accord with the remarks at the end of §4, t= Then (4. 14) 

and (4. 15) become 
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J wj 

= ^ L F^{t)cosu,^tdt 


0 

2ir/u) 


dB^ Zir/0,1 

■at^ = -rSS- I FJt)sinu,^tdt 


( 6 . 2 ) 


The linear coupling terms (4. 20) - (4. 23) all vanish when uu = no) , and the 

n 1 

only non-zero values of (4. 26) - (4. 29) are 


= 6 . 

1+ n, i+j 


4+ n, i+j 


The equations (4« 16) and (4, 17) are now 

^ 00 00 


=6 . . + 6 . . 
1- n, x-j n, j-i 

(6. 3) 

= 6 ... 6 . . 
4- n, i-j n, i+j 

(6.4) 


p=-yD A -y — B y [c ..a..8 . 

t 2nnn2o) n Z(jO Lt L mj ij n, i+j 


i=i j=i 


+ D .. b ..(6 . . + 6 .j^.) 

raj raj' n,i-j n,i+j'J 


(6.5) 


dB , 

-3^ = - i D B 
dt 2 nn 


” "i=lj=l 

- D .. d..(6 . . - 6 . ,)1 

raj ij' n,i-j n,J-l J 


( 6 . 6 ) 


With the mode shape >1^ = cosk t, the integral I .. is 

n n 


I •• = T (6 .^ +6 . . + 6 . .) 
raj 1 n, i+j n, i-j n, j-i' 

and the formula (3.60) gives eventually 


C .. = 
nij 


- ^ i+i 

4y n J 1 n, i+j 

+ — {[2a)^^-(y-l)(a),+w.)^] + (3-7 )(w.^-w.^)}(6„ . .+6„ . .) 

4y n J 1 J 1 n, i-j n, j-x 


(6.7 ) 


( 6 . 8 ) 



-34- 


o^. = - -4: {C 2 w^-(y-l)(a).- 0 ).)^J + (3-y)(co.^-co.^)}5 . . 

4y " J * J ^ 1^# 1+J 


+ {to^ + (w.^ - W.^)} (6 . . + 6 . .) 

4 y “ J 1 n,i-j 


(6.9) 


Let 


a = - tD 
n 2 nn 


e = - 5-^ E 

n Zco nn 
n 


( 6 . 10 ) 


and after some arithmetic. Equations (6. 5) and (6. 6) can be put in the form: 
•iA (3 * 

^ = o A + 0 B +~ > rA.(A .-A. -A _,.) - B.(B .+B. +B .)] 
dt n n n n 2 L i' n-i i-n n+i' i' n-i i-n n-i' 

i=l 

y .-(01^ -0}?)A. -(w\.-W.^)A .] 

2 L L,.2 1 ' n-i 1 n-i ' x-n i' i-n n+i i ' n+x 


• 1 ^ 
x=l n 


(6.11) 


- ^ B.[(W^ . -0).^)B .+(W. ^ -0).^)B. HU3^,.-(ji})Vj > 

^2 1 ' n-i 1 ' n-i ' i-n i ' i-n n+i i ' n+i J 

n 


dB 

1 

dt 


00 


n 


P V' 

a B -0 A +-^ ) [A.(B .+B. -B ^.) + B.(A .-A. +A ^.)] 

n n n n c u i' n-i i-n n+i i' n-i i-n n+i' 


i=l 


00 


+ -^ 1 A.[(co^ . -W.^)B .+(C0.^ -03.^)B. „-(<•) 4- 

2 L 2 1 n-i 1 n-i i-n i i-n n+i i n+i 


. , 00 
1=1 n 


( 6 . 12 ) 


+ -1 0 
00 

n 


where 



(6.13) 


The first series in (6.11)and (6.12) arise from the symmetric parts 


of C .. and D ... One can verify directly, in accord with the remark fol- 
nxj nxj 
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lowing (3. 6C) that, because terms cancel one another by pairs, the second 
series in (6.io)and (6. 11) vanish. Consequently, the equations to be solved 
for problems involving purely longitudinal modes are 


00 


® ® y CA.(A .-a. -a .+B. +B ,.)] (6.14) 

at n n n n 2 Zj i' n-x i-n n+i i' n-x i-n n+i' ' ^ 

i=l 

dB ^ 

= crB-9A y Ca.(B .+B. -B _)+B.(A .-A. +A ..)] (6.15) 

dt n n n n Z u i' n-i i-n n+i' x' n-i x-n n+i' ' ' 

i=l 

Some numerical examples are given in Section 10 and 11. For most cases 
considered, five modt's will be treated. The explicit equations, obtained 
from (6. 14) and (6. 15), are 
dA. 


T 




dB 


= ajAj+0jB.-p(AjA2+A2A3+A3A^+A^A5)-P(BjB2+B2B3+B3B^+B^B5) 

(6. 16a) 


1 ft 

■3T " ®iSr®iAj+2 - C(BjA2-AjB2)+(B2A3-A2B3)+(B3A^-B^A3) 

+(B^A3-B5A^)] (6. 16b) 


dA2 

"dT 


dBj 

"HT 


dA3 

~W 

dB. 


02 A^ + ^2 f 1 ^ 3 ^-^2 3-^5 ® f ® 1 ® s'*" ®Z ® 3 ^5 > 

(6. 17a) 

(6. 17b) 

a3A3+93B3+3 pAjA 2-3p(AjA^+A2A5).3p{Bi82+BjB^+B2Bg) (6. 18a) 


(6. 18b) 
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Z Z 

= a4A^+0^B^+4p(iA2 +AjA2-AjA5)-4P(iB2 +BjB3+BjB5> 
dB 


dAc 

-g^ = a5A3+05B5+5p(AjA4+A2A3)-5P(BjB4+B2B3 

dB^ 

-at- - «5B5-93A3+5 P(AjB3+BjA^+A2B3+B3A^) 


(6. 19a) 
(6. 19b) 

(6. 20a) 
(6.20b) 


VU. AN APPROXIMATION TO THE INFLUENCE OF 


TR/NSIENT SURFACE COMBUSTION 
Only the simpl;;st representation of the infl'ienco of unsteady com- 
bustion processes will be covered here. Elementary results for the linear 
response to harmonic pressure vaiiations form the basis. Certain of the 
features of truly transient behavior will be ignored in the interest of obtaining 
formulas which are clear and inexpensive to use. 

The influence f surface combustion is contained in R , defined by 
Equation (3. 42). It follows from (3. 40) that the corresponding contribution 


to the force F in (4. 1 ) is 
n 


,(c) 


J. ± 


n 


P E 
^o n 


2 at 


dS . 


n 


(7. 1) 


To simplify the discussion, consider only one of tlie pieces of ft ; e. g. , let 
6x = 1» fill = 0 » so (3. 42) is 

- at; 

' o 


ft 




(7.2) 


It is best here also to avoid the complications associated with a condensed 
phase: set k = 0 , so the following results apply to propellants not containing 


metal. 
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There exists a class of analysis, discussed by Culick (1968), which 
produces a formula for ih^/m^ , the fluctuation of mass flux due to a sinu- 
soidal variation of pressure. This is a linear result, m^/m^ being propor- 
tional to the pressure fluctuation: 


^ r\ 


(7.3) 


The response function, is a complex function of frequency, 

nAB 




X + y - (1+A) + AB 


(7.4) 


in which A is proport'.onal to the activation energy for the surface reaction, 
and B depends on both A and the heat released by the surface reaction. Out 

A 

of the same analyses, one can extract the formula for AT/T^ [see also Krier, 
et al. (1968)1 : 

^ = r M (R _ ILlil A (7 5) 

T L T C E b ^ y J p ' 
o o p ’ 

Consequently, if these results are used, (7. 2) becomes, for sinusoidal motions. 


o b p_ 


(7.6) 


Note that and are dimensionless functions of the frequency and the 


other parameters; if nonisentropic temperature fluctuations are ignored, R. 

rJ*'^(uu.) and r/^^ = R.j^^(t’.) - 
b 1 1 b 1 

Now the formula (7. 6) is, by construction, for steady sinusoidal vari- 


(r) 


ations only. The approximation suggested here is a means of using the formula 
under conditions when the amplitude and phase of the oscillations are varying 
in time. This is done by noting that for sinusoidal motions, i is equal to 
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u) ^ 9/ 9t . The replacement is made in (7. 6), and assumed to apply to all 
modes. Thus, p/p^ stands for ▼'jt- t and for an arbitrary pressure field ex- 
panded in the form (1. 6), R in (7. 1) will hereafter be taken as 




(7.7) 


The subscript ( )j on and means tiiat each function is evaluated at 

the frequency of the i^^ mode: these quantitie.-i can be calculated from (7. 6) and 
(7. 4). 

The force (7. 1 ' is now 



E 


n 


Qo 

I 

i=l 


<r.ds 

1 111 I « ^n 1 


(i) 


(7. 8) 


The further approximation, ii. « -m. n. , has been made in (7. 8): this is con- 
sistent with approximations already made in deriving the equations (4. 1). 


Finally, the rule (6. 8) gives directly the contributio~s fron» surface 
( c ) { c ) 

combustion to Gt an^i 0 to be used in (6. 12) and (6- 13): 
n n X / \ / 


a 


(c) 


n 





n 



a 


(c) 

n 


(7.9) 


(7. 10) 


1 

The same procedure c an be applied to propellants containing metal; only 
some details are chan.^ed to account for k ^ 0 . 

A similar approximation can be used for handling combustion within 
the volume of a chamber as one finds in liquid rockets, thrust augmentors, 
and rolid rockets exhibiting residual combustion. Although other contribu- 
tions will in general arise, associated with mass and momentum exchange, 
the direct contribution of energy release is represented by the terms con- 
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tainixg Q and in eq. (2.2). The perturbations are part of P* , (2. 11), 
which appears ultimately in H, and eqs. (3. 30) - (3. 32). Those 
functions are on the right hand sides of the oscillator equations (3. 39) and 
(3.41). The dynamical behavior of combustion witiiin the volume may, for 
example, be represented by some sort of response function; the well- 
known n-T model developed by Crocco and co-workers is a special form. 
In any case, die contributions to the individual harmonics can be approxi- 
mated as surface combustion was handled above. No results for bulk com- 
bustion have been obtained. 
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VIII. AN APPROXIMATION TO THE LINEAR AND NONLINEAR 
ATTENUATION OF WAVES BY GAS/PARTICLE INTERACTIONS 


Particularly ir solid propellant rockets using metallized propellants, 

but in other systems as well, some of the combustion products appear in 

the form of liquid or solid particles. The viscous interacting between the 

particles and the gas may, under suitable conditions, provide a significant 

dissipation of energy. It is often the case that the Reynolds number based 

on the particle diameter is outside the range in which Stokes' law is valid: 

it is necessary to use a more realistic representation of the drag force. 

This introduces another nonlinear influence in the general problem. 

Let dencte that part of F^rn eq. (4. 1), representing the influ- 

ences of inert particles. The terms involved are those containing and 

6Q' , the fluctuations of (2. 8) and (2. 9). By tracing the development from 
P 

(2. 10) and (2. 1 1 ) to (L. 18) and (2.19), to (3. 8) and (3. 9), to (3. 40), with H 
defined by (3. 30), one finds that the terms in question are 


(P) 


n 


5- f-J- [ ^ — 4r (6Q' +6u' • F H +6F' * ^ IdV . (8. 1) 

— at p p p^n p n-J 


p E 

n 


The differential heat transfer and force acting, per unit volume, between 
the condensed phase and the gas are defined by (2. 8) and (2. 9); explicit formu- 
las can be found only *>y solving the equations of motion (2. 5) and (2. 6; with 

the force F and heat transfer Q specified. Numerical calculations 
P P 

[Levine and Culick (197 2, 1974)] have shown that for many practical cases, 
nonlinear interactions are likely to be important. The approximate analysis 
here will be based on the nonlinear laws used in those works: 


-Ft — iSijt “•iTi 1 1 15 2/3 -I 


p 


( 8 . 2 ) 
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12C u 

Q. = p ^ 


2 (T -T)[l+ .23Pr-^^Re'®^] 


where 


P Prp^a^ P 


p 

Re = lu -u) . 

U P ' 


(8.3) 


(8.4) 


Hereafter, the real flow will be treated only in a local approximation 

so that the particle motions may be treated as one -dimensional: interactions 

between particles are assumed to be negligible. I'or a single particle, with 

spatial variations of tne motion ignored, the equations to be solved for u^ 

and T ' are 
P 




(8.5) 




dT' 


1 


12C Li I- *1 

P_2.(T' -T')Li + .Z JPr’^^Re J (8.6) 


p ^ Prp <7 

To evaluate 6F' and 6Q' , the quantities 6u' ~ u* -u* and 6T' = T’ -T* are 
P P P P P P 

required; by subtracting du'/dt from (8.5), and dTVdt from (8.6), one 
finds the equations 


with 


d6u' 

E 

dt 

.-L 

6u' = 

P 

du' „ 

■ dt ■ *^1 

6u' 

__E 

■^d 

2/3 

1%' 

(8.7) 

d6T' 

P 

dt 

. -i- 

6T' = 

P 

dT' 

dt 

6T' 
— E 

. 55 

Uu- 1 

F 

(8.8) 



”^d 

2 

P 

s 

18u 



(8.9) 




PrT 

2 C ^*^^d 
P 


(8. 10) 



^1 

t / P ^ 

- l( 

2/3 


(8. 11) 



^ \ u 

/ 
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K, = .459Pr‘^^( 


P ° 

-J£_ 

ut 


55 


( 8 . 12 ) 


Comparison of (8. 7) and (8. 8) with (2. 8) and (2. 9) (again with spatial varia- 
tions ignored) gives the formulas 



(8. 13) 


(8. 14) 


As an approximation, the nonlinear terms in (8. 13) and (8. 14) will be evalu- 
ated by using the line.ir solutions for 6u' and 6T' . With K, = K, = 0 , the 

P P It 

linear solutions to (8. 7) and (8. 8), satisfying the initial conditions 6u' = 6u' , 

P po 

6T' = 6T' at t = t . are 
p po o 

p t 1 

6u' = |_ — e Je ”u'(t' )dt'-a'(t)J -L6u' -u'(t^)]e ° (8.15) 

P d t P 

o 

6T' = -!-e fe ' T'(t')dt'-T'(t) -[6T' -T'(t)le ' (8.16) 

P U -I po O 

O 

The second parts, arising fronn the initial conditions, represent short 
term transients which are negligible for t-t^ these are retained, 

they will introduce in the oscillator equations (4, 1) terms which depend 
explicitly on the history of the motions. In the interests of simplifying the 
analysis, these term: will be ignored. This is an approximation which is 
accurate only if the p-^riods of the oscillations (al) harmonics) are long 
compared with t^. For the n^^ acoustic mode, the velocity and temperature 
fluctuations are 


u’ 




<^'Vn*n 




(8. 17) 
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where is given by (1. 12). Substitution into (8. 15) and (8. 16) gives 


1 

% = ^l^'^n - Vn^^:^ “HT 

n 


(8. 18) 


6T' (tl) T. Tl + '—) T i|t 

p (D'-ntn (1)/ — on 

n “ V 


(8. 19) 


The functions A and B are taken to be constant in this part of the calcu- 
n n 

lation because the results will eventually be used in the right hand sides of 
(4. 14) and (4. 15): becH.use the short term transients have been ignored, the 
lower limit on the integrals is t^ = 0. Explicit dependence on frequency is 


contained in Xj and X^- 


Xj = («^n^)/(nn^) 


Xj = («„nj)/(i+n, ) 


( 8 . 20 ) 
( 8 . 21 ) 


where 


a .Id t n t 


( 8 . 22 ) 

Substitution ol (8. 18) and (8. 19) into (8. 13) and (8. 14) gives for the 
linear parts only, 

di|r 


,5P.) . 

' p' — , L \ tu Tj n/ n uu nJ dz 


lin 


n 

p (V 

V 


"n d 


ri^n - _ 4 . ... 


(8.23) 


(8.24) 


n 


n t 


uu 


n 


These locally one-din.ensional results can be used in (8. 1) with di|f^/dz re- 
placed by Vi|,' to give 


(8.25) 
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Again by applying the rale (6. 8), one finds for the linear contributions from 


gas /particle interactions: 



(8.26) 

P 

? O 


n 2 1+H L ^ <1 J 

(8.27) 


Recent numerical results reported by Levine and Culick (1974) have shown 
that the result (8. 24) is quite "ood for smaller particles, and if the frequency 
is not too high. Beyond limits which are presently not well-defined, the 
Reynolds number (8. 4) becomes too large for the linear drag and heat trans- 
fer laws to be accurate. Further comments on the accuracy and some ex- 
amples are given below. 

The problem of analyzing nonlinear particle motions is avoided here. 

A correct treatment vould involve solving (8. 7) and (8. 8), the results then 

being used in (8. 13) tnd (8. 14). A very much sin. pier coarse is taken here; 

the linear solutions are used everywhere for 6u* and 6T' . The nonlinear 

P P 

part of the force F is 
^ n 


rF<»» 


] 

non! in 


p E 
^o n 






R/C r CK, , 




where the integrals are: 

2/3 

T = f Un- 1 6u' * Vllf dV 

1 p p n 

^4 ‘ 


(8.28) 


(8.29) 
(8. 30) 
(8. 31) 


(8. 32) 
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To evaluate these integrals, it is easier to re-write the formulas (8, 


and (8. 19) for three-dimensional i.iotions as 

V 1 JL 

6u' = ( > (-A +B ) sin{uj t+<j> ) Vft 

p \ — d / d n n ^ n 




n 


6T' = (- -)(l+n.^)^(A^+B^)^cos(u) t+4)‘^)t|f 

p yO'Ct'/ t 


n 


n n ^n 


where 

sm 4> = 

n 


(B^-n ,A^) 

n an 

? 1 2 7 7 * cos (j) 

(i+n/)2(A +B )2 
d n n 


(A„+n B ) 
n cl 




sin 4> 


. (B -CI A) 
t n t 


' t n n 


. (A„+n, ) 

^ . t n t n 

cos (p = T* 1 — X — ■ * ■ ' — 


The formulas (8.29) - (8. 32) become 

1+1 


1+e, 


* X ^1 ^ bj 

I = I k (l+n,^) ^ ( — ^-) (A^+B^ ) ^ sin((i! t+<t>^)jsin((jj t+<j>'^)| 

llnci \— / nn nn nn 

yK 


X , \ 2 / — 


y\ V n 


1+^: 


d 


{-cos(uu t+c])^) 1 sin(tt' t+<(>^)l 
^ n n ‘ n n ' 


3 3L~- J n ndt^- n n 

y k 

9^L 2+f;, 2+e, 


1 

2,2. 


r 


(i+n, )‘"x 


^4 = '4L— r 


Vk 


i”| (A^ + B^) ^ {sin^(o} t+r I sin((j) t+<i)^) I ] 

J n n dt^ n n’ n n J 


n 


where ^ t = 55 and the integrals involving the mode shapes 


J(V4_^)^|vtJ >dv 


n 


ly = — ^ Jc ^ Ivt 1 ^dV 
2 *J n ' n 


18) 

(8. 33) 
(8. 34) 

(8.35) 

(8. 36) 

(8. 37) 
(8. 38) 

(8. 39) 

are 
(8. 40) 


(8.41) 
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-T J 


(8.42) 


‘ "dV 


(8.43) 


For comparisr-i with numerical results, these integrals will be evalu- 
ated for longitudinal modes in a uniform chamber. Then dV = S^dz and 

A = cosk z with k mr/L and 0 ^z S L. By making use of the periodicity 
n n n / o t- 

of ilf . and with 0 = k z , one c jn show that the values are 
’^n n 


3 +e,. 2 


n 0 


2nS tt/ 2 2+1, nS 2+|, [r(-^')] ,nS . 

= ‘ nUT : " - -Kt^) 


I = -j- J cos 9sm 0d9 = -g 

„ 0 n r(_^) 

S mr 2 

I = J cos 0 sin 0d8 = 0 


,nS 

..644(^) 


(8.45) 


(8. 46) 


S nv 


am M I 

I = J COS 0 I sin 9 1 ^d0 = 0 


(8.47) 


To simplify writing, assume only here that t^. = which is very closely true 

d t 

in many practical casos, so cf)^ ~ ^n* define the functions 


C (t) - sin( n t+4> ) sin(ou t+4> ) 
n n ’ n n ' 

Kn<‘> = / 


(8.48) 


(8. 49) 


With all the preceding brought together, the nonlinear part of the 


force (8, 1 ) is 


[F^P^ =ui) (f~){w ,(l+fij)^A^ + B^7^u) C (t) 

i; n 1+K ^ nl d n n n n 

nonlin - 


(8. 50) 



rfpkoducibiiOT ot'M 

'» 1NA.L PAG® ® POG“ 


where 


C VO) . n t 



■^1 'I 

y > 

^ ID Cl ,/ 

Y«>„ 

n d 

.1)~K 


C 

P 



(8.51) 


(8.52) 


The time-averaging of (8. 49) according to (4. 14) and (4. 15), with 
T= Tj (bee the remarks at the end of §4) requires the integrals 

2tt/o), cosuu t !>_ fl 

, 1 , n 1 , 2u'n r cos0 


5 Vl ^ 


H / UUl V.V/OUU Ir 

^ 1 ^ n ^ j ZTTn j CO 

T, '^n^n^ -sinu) U ~ Zir ^n"i -si 


2ix/aj. 


COS(D t 


6 'L' 


J_ J * J . "Jd, = ± 

T, n^i -sinij) tJ 2 tt 


sin0 / 


Zirn r cos0 


n 1 0 


[ ^nl.sinel^® 


where 9 = J-'^t . With and given by (8.48) and (8.49), the integrals are 


sin(0 + 4.^)| ^ d0 

^ ' { !sLV}^ 


Because the .tegranrts have period 2 tt , the limits on the integrals have been 
changed from (0, 2Trn) to (0, 2tt); the integrals are then multiplied by n, giving 


the factor n/fi t, - (^.tt 
n 1 


Now let il/ - 6^ 4^^; the limits become (4>, 2rr+4)) 


which, again because of periodicity, can be changed to (0,2ir). Moreover, 
it is easy to show that the integrals over (tt, 2'rr) are equal to those over (0,tt), 
so one has 

- TT ^cos'i’coscb +sini'sin<J) 1+4^ 

r I r ’ n t . , , , 

M • ' J i ^ I - d’!' 

5 ^ X -sinvccsc), -f- cos vs in 4) j 



-48- 


I 


^ jT.ros^c.84.^+sin^rsin,i,^ , , 1 + ^ _ 

6 "it Ij I -sin\|icos(j)^+cosil;sin4)^J t ^2 


e,-l 




All integrals containing cos(l; vanish, and those containing sinil/ can be re- 
duced further to the interval (0, tt/ 2), giv ng the result?. 


, , sin4. , W2 . 2+e, 

's = 1 -cos* ; I ““ ♦ 

n ^ 


d^lf 


y psincf) > p TT/^ ^9 “I 

h " Fl-cos'i 

n 0 0 


-t/2 2+1, V2 I: 

0 “0 
For = 2/3, the integied. in has the value 0. 698 as in (8.44), auid for 
^2 = 0- 55, the integrals in are 


tt/ 2 2-h|2 

f sin il'dili =0,713 , 

•I) 

Finally, then, the inti^grals are 


V^' "2 

1 sin ilidii' = 1- 175 

'0 


.. 396 r 1 

IT 1 -cosd) / 


(o. 53) 


I = > 918 r 
6 TT L -cos (j) 


( 8 . 54; 


After sin())^ and ccs4>^ are replaced by thrir definitions (8. 35), one 

finds for the contributions to A and B , from nonlinear gas /particle 

n n ^ 


interactions: 


dA 

~ar 

dB 


r = ( 8 . 55 ) 


= IK (-r^)^V ,(-A -fi,B )?V ,(-A -n,B )| (8.56) 

dt n\ l + K/ ^ nl n d n nZ n t n J 


where 


V . — 2 ^ - w ,( A ^+ B ^ 

n 1 TT nl n n 


.353 


'8. 57a) 





TT nl n n 


Note that in (ft. 55) and (8. 56), appears explicitly in two places, and in 
two places. Actually, what arises is 0= x^t where, because of the assump- 
tion preceding (8. 48), t = ^ arbitrary choice t = or t r 

has ’ oen made whi suitable to simplify somewhat (8. 55) and (8. 56). The 
numerical errors incurred should normally be very small. 

A few results have been obtained for the attenuation of a standing 
wave initially excited in a box of length L. containing a gas /particle nriixture. 
These have been carried out for direct comparison with the numerical 
results reported by Levine and Culick (1974). First, two cases wilx be 
given in some detail. The particle diameter is 2. 5 microns, the particle 
loading is = 0. 36 and the frequency based on the equilibrium speed of 
sound a is 800 Hz in each case. The material and thermodynamic proper- 
ties are lii>ted below: 


Specific heat of the gas 

Specific heat jf the particle material 

Isentropic exponent of the gas 


Prandtl number 


T emperature 


= 2021,8 Joule/kgm- K 

C = 0. 68 C 
s p 


y =1.23 
Fr = 0. 8 
T = 3416 


Pressure 


/iscosity of the gas 


= 500 psi 


M = 8. 834x10"^ kgm/i 


Density of the condensed material 


= 1766 kgm/m' 


The only difference between the two cases is that for one, the initial 
aniplitude is Ap/p^ = 0, 03 and for the other Ap/p^ = 0. 15, 
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Figures 8-1 and 8-2 show the waveforms obtained from the 

approximate and numerical analyses. The approximate results are 

the solutions to Equations (6. 16)-(6. 2C); the linear coefficients a , B 

n n 

are given by (8a 26) and (8. 27) and the additional nonlinear terms by 
(8. 55) and (8* 56). The period used to normalize the time scale is the 
period 2L/a for linear waves based on the speed of sound of the mixture. 

The similarities betv/een the waveforms and the generation and decay of 
the even harmonics is apparent. There is a difference in the relative 
phases between the even harmonics and the total waveform. This seems 
to arise almost entirely in the first cycle of the oscillation; it may be 
due to details of the numerical routines and the way in which the compu- 
tations begin. This difference may therefore not reflect a genuine differ- 
ence between the approximate and “exact** analyses. 

A more quantitative measure of the behavior is the instantaneous 
**value of the decay constant, a^, calculated for successive peaks of the 
waveform. The histories of for the two cases are shown in Figures 
8-3 and 8-4, Further remarks on the behavior of or^ may be found in § 7 
of the report by Levine and Culick (1974), The purpose Vi3ie is only to 
demonstrate that the approximate analysis does give fairly reasonable 
results for this case. There are, however, limitations, not yet clearly 
defined, which arise from the approximate treatment of the nonlinear 
acoustics as well as the gas /particle interactions. 

The linear behavior used here, described essentially by the velocity 
and temperature lagF given by (8. 18) and (8. 19), is, in fact, quite restrictive. 
Although it is not ap].arent from the analyses pre sented here, the work 
reported by Levine and Culick shows that the results (8, 18) and (8, 19) 
are accurate only when co is small. If this condition is met, then 6 Up 
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and 5T^ are relatively small, which is consistent with the general nature 

of the perturbation analysis used here* For the 5-micron particles. 

jiT^ ^ * OS SOO Hz. The waveform and Cl for 10-micron particles at 
d P 

1500 Hz ^ 2.4) are shown in Figs. 8-5 and 8-6 for an initial ampli- 
tude of 3 per cent, computed with the approximate analysis. For this case. 

^ -530 sec*^ when the amplitude has decreased to 1 per cent; the more 
exact numerical calculations give ci^ «870 sec*^. Finally, in Fig. 8-7, 
the waveform calculated with the approximate analysis is shown for 30- 
micron particles at 1500 Hz « 25. 8). The vaxue of is -70 sec ^ at 
an amplitude of 1 per cent, and the numerical analysis gives -148 sec 
Note that the fractional error between the approximate and numerical re- 
sults for dp increases with increasing (ut^. The greater amount of har- 
monic content present when larger particles are considered is due mostly 
to the reduced drag and hence reduced attenuation at the higher frequencies. 
The results from the approximate analysis can be improved for highe:*: val- 
ues of by using different functions X^.X^ instead of (8.20) and (&. 21). 

Wave propagation in a gas /particle mixture is dispersive; the speed 


of propagation depends on the frequency and particle size, through the pa- 
rameter OL T^, as well as on the amount of condensed material present. 

For the problems treated here, the length of the box. and hence the wave- 
length of the waves, is fixed. Consequently, a change in the speed of sound 
it. reflected as a shifr of frequency. The periods of the waveforms shown 


in Figs. 8-1, 8-2, 8-5, and 8-7 are not exactly equal to the period based on 
the equilibrium speed of sound. It is obvious from the figures that the fre- 
quency shift increases wit! In §12. it is shown that to first order, the 

frequency shift due to linear dispersion is 6f = -6/2ir ; here, B is given by 
th 

(8.27) for the n mode. The actual frequency change is dominated by 0^ ; 
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TABLE 8-1 


Summary of Results for Waves Attenuated in a Gas/Particle Mixture 


Particle 

Diameter 

(Microns) 

(Hz) 

'f 

(Hz) 

f 


6 f 
(Hz) 

.li. 

2ir 

^"p^Un 
(sec ^) 

(sec ) 

2.5 









Approximate 

800 

954 

800 

00 

o 

• 

1 

.74 

-59 

-61 

Numerical 

800 

954 

801 

00 

o 

• 

1 

— 

-58 

-59 

Aoproximate 

1500 

1788 

1505 

. 15 

5 

4.8 

^203 

-204 

Numerical 

1500 

1788 

1504 

. 15 

4 

— 

-201 

-204 

10 









Approximate 

1500 

1788 

1695 

2.4 

195 

189.9 

-509 

-530 

Numerical 

1500 

1788 

1740 

2.8 

240 

— 

-509 

-870 

30 









Approximate 

1500 

1788 

1740 

21.6 

240 

224. 1 

-67 

-70 

Numerical 

1 

1500 

1788 

1787 

25.8 

287 

— 

-69 

-148 


t 

f 

6f = 


^"p^lin 


frequency based on the equilibrium 
frequency based on the speed for the gas only, 
frequency of the calculated wave 
f - f 

e 

decay constant for linear waves, Eq* (8.76) 


a 


g 


[(y-i)CpTi^ 



decay constant for the calculated wave at approximately 1^ amplitude 


NOTE : Small differences arise in some quantities (CU T ^ and ) 

which should have the same values when calculated by the approximate and 
numerical methods. This is due to a small difference in the value of some 
property, most probably the viscosity. 
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TIME/PERIOD 

Fig. 8-l{a) Attenuation by 2. 5-micron particles at 800 Hz 
according to the approximate analysis, Ap(0)/p 
0. 03. ° 
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Fig. 8- 1(b) Attenuation by 2. 5-micron particles at 800 Hz 

acc( rding to the approximate ana lysis, Ap(0)/p = 

0. 15. ° 
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time/per 1 00 


Fig. 8-5. Attenuation by 10-micron particles at 1500 Hz 
according to the approximate analysis, Ap(0)/p 
= 0.03. ° 
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Fig. 8-6. The decay constant for the case shown in Figure 8-5. 
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Fig. 8-7. Attem ation by 30-micron particles at 1500 Hz ac- 
cording to the approximate analysis; Ap(0)/p^ - 0. 03. 
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values of and 5x calculated from the waveforms are given in Tab'e 
8-1, a summary of the computations discussed here. 

These, results serve to demonstrate that the nonlinear generation of 
harmonics has substantied influence of the detailed character of the attenua- 
tion of waves. The rather close agreement between the approximate rasu’ts 
and the numerical calculations under conditions when the approximation to 
the gas /particle interaction is more accurate suggests that the approximate 
treatment of the fluid mechanics (i. e. , the terms represented by I in (3.45)) 
is realistic, at least for moderate amplitudes. 

Fi om the point of view of reducing data for the attenuation by 
particle/gas interactions, it is annoying that the veJue of the decay constant 
changes so much as the waves die out (cf. Figs. 8-3, 8-4 and 8-6). It is 
possible that if the higher harmonics are filtered from the waveform, the 
behavior of a for the first harmonic alone may not be so extreme. Calcu- 
lations have not been done to check this point. 

Perhaps the Simplest check of the fluid mechanics alone is calcula- 
tion of the behavior of a wave in a box with no particles present. In this 
case, the wave must of course steepen, eventually forming a shock wave. 
Neither the exact nor app’"oximate analyses can accommodate strong shock 
waves, but the initial period of development may usefully be examined. 
Figure 8-8 shows the exact result, and Figs. 8-9 and 8-10 show the results 
of the approximate analysis when five and ten modes are accounted fo. . 
Again, the qualitative agreement is quite good. As one would anticipate, 
the period of the wave decreases as the amplitude increases. For both the 
approximate and numerical analyses, obvious distortion of the peak o.curs 
at about the fourth cycle. (The sharp jags in Fig. 8-8 may be due in pnrt to 


the numerical routine. ) 
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Fig. 8-9. Steepening of a standing wave in a pure gas according 
to the approximate analysis, with five modes accounted 
for; f = 900 Hz. 
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Fig. 8-10. Sti epening of a standing wave in a pure gas ac- 

coidinp to the approximate analysis with ten modes 
ac< ounted for; f 900 JIz. 
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IX. AN APPROXIMATION TO NONLINEAR VISCOUS 
LOSSES ON AN INERT SURFACE 
There are two reasons for examining the influence of viscous 
stresses and heat transfer at an inert surface: these processes can be 
significant stabilizing influences; and in the presence of oscillations^ 
the average heat transfer increases substantially* The second is a 
nonlinear effect wh:ch has on occasion caused serious structural prcb- 
lems» particularly in liquid rocket motors* The main purpose of this 
section is to show one way of incorporating some experimental results 
within the analysis developed earlier* As part of the argument, the more 
familiar linear resi’lts will also be recovered. 

The viscous stresses and heat transfer at a surface are, of course, 
associated with a bc^undary layer, but they can te accommodated here by 
suitable interpretation of the force F and heat source Q in the equations 
developed in § 2 , It is only those terms which are required in this 
discussion, so the wave equation for the pressure fluctuation is simply 



-2 

a 



R 3Q' 
C at 




The boundary condi Ton is 


( 9 * 1 ) 


n* Vp* = n* F* 


( 9 . 2 ) 


For this problem, then, the equation for the amplitude of the n^^ harmonic 
is 

o a E G 
^o n V 

The heat source Q' is taken here to be associated with tha heat flux vector 

** 

q*. and the force F with the viscous stress t-^sor T: 
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Q' = V.q* 


F* = 7- T' 


(9.4) 


Both q' and T' are Fignificantly non-zero only in thin regions near the 
boundary. Then if y denotes the coordinate normal to the wall, measured 
positive inward, dv = dydS and 


Q' - 


9y 


0T • 


(9. 5) 


The conventions used here are that q* , the component of q* normal to the 
surface, is positive for heat transfer to the wall, and F*, being parallel 
to the surface, is positive when the force tends to accelerate the gas* In 
the volume integrals of (9. 3), and are essentially independent of 


y, so one can write 

00 9q * 

/ Q’ dv / / dS >1^ / dy = - // 


n 


n dy 
o ^ 


'C+n'*® 


_ 00 0T ' 

/ F'. V i dv J I dS Vijr • / dy = -If T ' . V\|f dS 

n n 0y ^ w " 

o ^ 


n 


(9.6) 

(9.7) 


Here, the surface stress is = - M(^’/9y)^i where u* is the velo< ity 


fluctuation parallel :o the surface, so 


-jfT'-Vf dS = //(/l¥-) dS 

w ^n ' 9y 'w ’^n 


(9.8) 


Only the linear stress will be treated, but both linear and 
nonlinear contributions to the heat transfer will be accommodated by 
writing 


3 ' = (k ^) + h(T ' - T ' ) ;9.9) 

^ dy w 00 w 

Owing to the thermal inertia of the wall, » U; the temperature fluj* 

tuation T * far from the wall is that associated with the acoustic field. 

00 

The average heat transfer coefficient, h, will be assumed in a manner 
described below, to depend on the amplitude of the acoustic field. 
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equation (9# 3) can now be written 


T + OJ T1 ^ 

n n n 

p 

n 


"o n 


P E** 
*^o n 


V 

,2 dt 


// h T ' ^ dS «9. 10) 
00 n 


First, calculation of the linear parts will be outlined. The known 
solutions for the velocity and temperature fluctuations in a sinusoidal acoustic 
field are 


u' = 

where u, T are the amplitudes far from the wall, and 


l-d+i) 


6 = J 2v/fa* 


Thus, for purely harrronic oscillations. 


(9. 11) 


(9.12) 


/ , 
«*8r’w 

9T* 

' 9y ' vv 


u 

|(Pr )5(I«)T e*“‘ 


As in the preceding section, the replacement i co 3/9t is made, 
th 

and for the n harmonic one has 


au* ^ 


i (:~2 ) "n ^ I 

o \ y ^ c / n n n nil 


(-W-) = 


where (V ) .. is the component of V parallel to u* at the surface, Sub- 
n ” ^ 

stitution into (9. 10), and some rearrangement leads to 
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TJ /Q 

’n"“n ’V' V- ^ 

p iL, y 

n ' 


(W„v/2)2 


(V) . '“'-'"^ 1 _L //r/!!!L\' t // *2 

n 2(l+c ) „2 ■'•'LVk / ^ V 


ITj' E 


-/Pr 


• 9.13) 


«9. 14) 


it has been assumed that the motion far from the wall is isentropic, so 
r^-(y-l)(T^p*/ yp^ )• For a longitudinal mode in a straight cylindrical 
tube, 

V y 2 

//(-^-) dS = // dS = irDL/2 

and with C =0, (9. 1 1) becomes the familiar result for the decay con^tant 
m 

:or a standing longitudinal wave: 



The treatment of the remaining term in (9* 13) rests on appeal to some 
recent experimental results. Perry and Culick (P.74) have reported meas- 
urements of the time- and space -averaged values of the heat transfer c >ef- 
ficient in a T -burner with propellant discs at the ends. Denote this coef- 
ficient by (h ). The data could be quite well represented by the expression 


= 0.044 Re j (9.15) 

K a 

where 

IpI 6 

Re^ = ^ (9.16) 

H a 

The symbol lp| ^ represents the maximum amplitude of the pressure iluc- 

tuation, namely that measured at the end of the T -burner; this is equal to 

p It) I for the n^^ mode, 
o n 

In the absence of any other information, two assumptions will be 


made. First, it is reasonable to assume that the local time-averaged heat 
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.ransfe** coefficient is also proportional to the square root of the pressure 

implitude, as (9. 15) shows for the space -averaged value. Hence, for the 
th 

a mode. 


h = K U Ti 1^ 
n ' n n * 


(9.17) 


where is a constant to be determined. The second assumption is that 
■ 9. 17) is valid for all modes. From the definitions of h and <h ) , 


<h) = i- //hdS = ^ / hdj 
^ ^ o 


(9.18) 


for a cylinder. It follows from (9. 15)-(9. 18) that the constant is given as 

lo/p / w \4ri L I t1 
K = 0.044 — f 
" TiTf '' Zv' l-L o ' -1 


The integral has the value 

r / K h dz 
1 j o ® 


1 L 1 

f / |cos k zj2 dz - .765 


,(0 \ir 

= 0.C575 2. (JL) 

JJTk ^ 2 v ^ 


(9.19) 


Note that the mode shape cos (k^z) is used to obtain (9. 19) because the data 
was taken in a uniform tube. The assumptions introduced above imply that the 
result is supposed to be valid for a local surface element whatever may be 
the mode in the chamber. 


With (9. 17), th:> last term of (9. 13) is 
^ \ T K . . .1 




The constants can be combined as 


H = 0. 0575 TT 
n 


.4 _y..iL 

Pr y 


(V f a ^)^ J 
n n 


( 9 . 20 ) 
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n 




n ' n 


(9.21) 


n 


Equation (9. 13) is now 


T) T] = -CO T] ) - H (ri Iri J^) 

^ ' n dt n n * 


n n n 


n n n n 


( 9 . 22 ) 


When substituted into the formulas (6. 1) and (6. 2) obtained with the 


method of averaging, the linear terms give 
2-rr/a:, 


, 1 , cos CO t ^ ^ A -B ^ 

^ { -sin nTt 1 <*' = 2 { a" } 

o n n n 


(9.23) 


Define the angle as 


B 


sincp 


n 


n 


V n n 


cos cp 


n 


n 


^ n n 


f9. 24) 


and one can write 


3/4 


^ h h = (A'^+b'^) sin(co +cp )lsin(co t4<p )\^ 
n* n’ ' r n n n ' n n 

The integrals (6. 1) and (6. 2) arising from application of the method of 

averaging to the nonliiiear term are 
ZttCO 


4- w h p)dt 

Tm _ I sinco^tJ dt n n' 

2 Tm r c j i 

i { -sine) 5e 


Zto 

n n ' 


Zim 


where 0 □ CO^t. The integrals can be reduced in much the same way as those in 


§ 8 to give 


/,2^^2.3/4 it/2 

(A. +B ) r coscp . ' 

^ ^ j sinca^ f / cos ^ sin^>|^d\^ 

^ o 


Ztt 


The integral has value . 478 and the final result for the nonlinear term in 
(9.22) is 
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2ir/co, , . TA (A^+B^)*^ 

^ ‘ «n * -smut} O' = B (A^B^)* r 

O n ^ n' n n ' ^ 


<9.25) 


The equations for and B^, with only the viscous losses shown are 


therefore 


dA 

n 

- - ■ 

dt 


dB 

n 

- - ( 

dt 

- % ' 


n n 


I A (A^ + )* 

n n n n ' 


(A +B ) - .457 H B (A^+B„^)^ 
' n n n' n n ' 


(9.26) 

(9.27) 


In the special circumstance when these equations are applied to waves in a 
cylindrical chamber w^th an inert lateral boundary. 


^ _ 1.836 

‘^n R 

c 

where is the radius of the chamber. Then the coefficient in (9. 26) and 
(9. 27) can be written 


.457 H = 
n 

where the units are: 


1. 143 y-1 

Pr y 2 

R , meters; 
c 


jq-4 1000 ' ^ 

_ 2 
V, (meters) /sec 


a .2^ 

Tm' J 

; f , Hz; and a , 
n 


(9.28) 


meters /sec 
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X, CHANGE OF AVERAGE PRESSURE ASSOCIATED 
WITH UNSTEADY WAVE MOTIONS 


There are often circumstances, particularly in solid propellant 
*Aiotors, when a substantial increase of the mean pressure may accompany 
unstable oscillations. This can become a serious practical problem, but 
it will be treated here only to the extent that one aspect is related to the 
analysis developed in the preceding sections. It should be remarked in 
passing that the likely cause of large DC shifts of chamber pressure in 
motors is usually a change in the ourning rate of the propellant, perhaps 
associated with the erosive action of unsteady motions parallel to the burn- 
ing surface. That subject will not be discussed h< re. 

The main point of this section is that the zeroth term in the expan- 
sion (1.6) represents a DC shift of pressure; the corresponding frequency 
and mode shape are = 0 , = 1 . In all the preceding discussion, and 

in other works as well, the influence of a small a' erage change of pressure 
has been ignored. All the formalism developed in sections 2 and 3 is valid 
and can be used to compute the change of average pressure due to the un- 
steady motions, as well as its influence on the wave motions. From (4. 1) 


and (4. 2 ), the equation for is found to be 



oo oo 


S (D .n.)- S D [A ..fi.r.+B ..'l.T.l 

i=0 " UO j=0 " J oij 1 i 


(IJ. 1) 


The values of the coefficients E^^ depend on the processes considered. 

For example, with D . given by (3. 41), one can easily determine the results 

111 

— , w 

D = ^ dv (10.2) 
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— , w 

D . = fi|(. -£■ dv 

Ol V *' 1 — 


(10. 3) 


no 


1 r H — 

\ J U • V i(r^dv+(v - 1 ) J* i|(^ dv 


(10.4) 


n 


where the factor v ar-ses from = v . It is interesting also to examine 
the term arising from surface combustion; this can also be interpreted to 
represent other sorts of processes as well. The term in question on ti e 
right hand side of (3. 40) is 

JJr ('^dS . (10.5) 


V 9 

p ^ 


For purposes of illustration, the fluctuation of the gases leaving the sur- 
face is assumed to be simply related to the pressure, and R is approxi- 
mated by (7.7 ), 

a. 

r. 

i=0 -1 

Equations (10. 5 ) and (10.6) give 
oo 


a 


, u, f: . 

ob^.nl-l UJ. 1 OtJ ll 


( 10 . 6 ) 


7 JJ' Tf . -uu.R.^^^n. 1 dS 

' b 1=0 1 'ill 1 ' '^i^n 


(i). 


(10.7) 


in which the approximation « -ou^ T|^ has again been used. Comparison 
with (3.40) shows that the corresponding additional terms on the right hand 
side of (10. 1 ) are 


/ 1=0 1 1 11 11 


1% 

V 


/ \ V U, oo .V ... 

ff R ^’^\i!;=0)ti dS 4 — JJ . (10.8) 

‘J'-'o O 1 111 11 


The first term represents the contribution to the change of average pr«»ssure 
in the chamber (really the second time derivative) because the surface 





combustion itself responds to the average change of pressure. 

Suppose, for example, that all other terms are ignored and that 

(r ) 

the response is ^iniform over the surface. Ihen (10.7) and (10. 1 i 

give 


^2 


dt 


y ^ 


b Q (r), 
— fto 


dn 

■aF" 


(10. 9) 


where is the total T.rea of burning surface. According to the discussion 

( r ) 

in ^7, see eq. (7.4) and following remarks, R may be approximated by 


the value of R. 


(r) 


at ' 
d^r 


dt 


= 0 ; thus, (fi;:^0) n , and (10. 9) is 

■ = ~r- "-ST ■ 


The first integral is 


dp’ 

^ o 

^dT 


- - "b . 

y n u, p' 

’ b V o 


(10. 10) 


where p^ is tlie change of average pressure; the constant of integration has 
been set equal to zero which merely sets an acceptable initial condition. 
Now (10. 10) is simply an expression of the change of pressure in a closed 
chamber due to mass addition. To see this, note that conservation of mass 
provides 


dt ^ 


where is the dens-'ty of solid and r is the linear burning rate. The 
perturbation of this e<'uation combined with the isi ntropic relation p 
gives 


p V dp' 
o 

”p dr“ 



But r' n(p'/p^^)r and p^r Isist equation becomes (10, iO). 

The point of this elementary example is lliat the elaborate formadism do- 



-76- 


veloped for unsteady motions does indeed contain familiar quasi-steady be- 
havior. A term representing the influence of erosion can also be incorpo- 
rated in this analysis, but it will not be considered here. 

It may be noted that the rema: aing terms in (10. 7) sometimes con- 
tribute nothing. For .!xe case of longitudinal modes with a burning surface 
extending the entire length of the chamber, for example, the integrals of 

\L'. over the surface vajiish. Interactions between the wave motions and 
1 

surface combustion th^n contribute nothing to the average pressure in the 
chamber. 


The nonlinear terms contained in 1^ , eq. (3.45), also contribute 
to the average pressure in the chamber. Only the parts containing p’ as a 
factor are involved, and it is not difficult to use the results (3. 58) and (3. 59) 
to determine the following formulas for the coefficients: 


B 


on 



B 


nen 


y - I 2. 

~ JL' 

— n 

V 


( 10 . 11 ) 


( 10 . 12 ) 


All others vanish; there are no non-zero values o' A .. for n, i, or j equal 

nij 

to zero. According U (10. 1), the nonlinear interactions produce a negative 
contribu" n to the se< ond time derivative of the average pressure: 


( 10 . 15 ) 


/ dV 



^ v-l ' 

V oo . / E . \ 

f ) 


. ( 




V - . 

/ iT 1 i V V / 


It is perhaps surprising that the sign is negative. Because in some sense 
the processes involved represent a dissipation of energy from the wav' mo- 
tions to the average s' ate, one might have anticipated a positive sign c .w- . 
sponding to a tendency to increase the pressure. The interpretation ij- not, 
however, clear at this time, 

Thi' influence of changing average pressure will not be c onsidi red 


further in this rcpori. 
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HCTHODOCIBILnY OF THB 
PAQ« B POOR 


XL PPPLICATIOX TO THE STABILITY OF 


LONGITUDILAL MODES IN MOTORS AND T-BURNERS 
It has been shown in §6 that the nonlinear terms simplify consi^ier- 
ably for the case of longitudinal m'^des. This situation will be treated 
here for several examples. Some preT ninary results were reported >y 
Levine and C’xlick (19!^4). Application to large m< tors was examined by 
Culick and Kumar (1974). The discussion here 5s to demonstrate h -w the 
approxinfiatc analysis can be used to study pr^ ticcJ configuratiens anc to 
provide a limited comparison with the more exact numerical results re- 
ported els e’.v here. 

For all the ca.culatioi s orted in this section, the following na- 
terial and thermodynamic pi '^-perties are used: 
specific heat of the gas 


c = 2021.8 
p kgm-^K 


speci- : heat of the p« rti< ’ material 
thermal diffus»vity of the propellant 
Prandtl number of the gas 
viscosity of the gas 
isentropic exponent o: the gas 
linear barni g rate oi the propellant 
density of the propeil mt 
density of the c^nden. ed material 


C ^ 0. 68 C 
s p 

K - 10“^ m^/sec^ 


Pr 0. 8 

- 4 , 


. 66 


u- 8. 834X10'^(T/3485) 

m-sec 


Y ^ 1.23 


..3 


r ^ 0 00812(p^/50C)*' m/sec 
z: 4000 kgm/m^ 
p = 1766 kgm/m"* 


^ ^ ^ Applica^ Lon to a Small Cylindrical Motor 

Because the n mlinear acoustics terms represented by , eq.(3.45), 


arc given explicitly, .he first s^ep in tha analysis consists in evaluating the 

contributions to the linear f^icients Ct end 9 , Four contributions to 

n n 

linear stability are ii cluHr J; arising f^-c m the nozzle, the condensed i 'la- 
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terial in the gas phase, the surface combustion pr< cesses, and the one 

dimensional approxim tion to inelastic acceleratioi of flow issuing from 

the lateral surface ("fJow turning”). The formula.' for the correspond! ig 

values of a are 
n 

Nozzle 


111 



Particles 

a 

n 

Flow Turning 




.l=X 


n 



( 11 . 2 ) 


(11.3) 


Combustion 


i V u 


b - V / 


(r) 


(11.4) 


The formula (11. 1 ) is based on quasi-steady behavior; the velocity fluctua- 
tions are in-phase wiiii pressure fluctuations at the nozzle entrance. Hence, 
the response function has no imaginary part and, as shown by eq. (7. U ), 

the value of for th * nozzle is zero. The term represen^^ing flow-ti.rning 
n 

is the last one in (3.41); there is no corresponding , so 9^ is als<, zero 
for flow turning. Tiie only non-zero values are for the gas /particle incer- 
actions and combustic n; the first is given by (8. 2' ) and the second by use of 
(7. 10) and (11. 4); 

Particles 


9 

n 


1 

n 

2 


\r 

K \\ d 




- r 


n 


(11.5) 
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Combustion 




( 11 . 6 ) 


The two examples treated here were discussed in ?7.4 of Culic c and 
Levine (1974), where che result of the numerical i_nalysis is given. E^.ch 
is for a motor having a cylir irical bore and lengtl. 23. 5 inches. The i lean 
pressure, port area, and th>*oat area differ and aie given below: 

(a) (b) 


Length (in. ) 

23. 5 

23. 5 

2 

Port Area (in ) 

3. 33 

4.73 

2 

Throat Area (in ) 

. 439 

.562 

Pressure (psia) 

1568 

1412 


These are the first and last cases given in Table 7-3 by Culick and Levine 
(1974). The fundaunental frequency is 900 Hz. 

For both cases, the particle diameter is assumed to be 2.0 microns 
and the mass fraction is k - 0, 36. The combustion response is taken to 
be the representation (7.4) given above, with A = 6. 0 and B - 0.56. be- 
cause the pressur :s ore differeiit, so are the flame temperatures in tr.e 
two cases; for (a), T = 3525°K and for (b), = i515°K. The small «if- 

ference has only minor influence on the results. 

With the abov^' values, and the formulas (il.l) - (11.4), waves 
for case (a) are founo to be stable. The values of the decay constants for 
the first five modes i.re given in Table 11-1 belo\ , The numerical Cc.Lcula- 
tions produced an unstable wave which, with an initial amplitude of 5 per 
cent (fundamental me le) eventually reached a limiting ampUtude of 4. I per 
cent. In view of the su^^essful comparison for the cases treated in § 8 for 
attenuation by particl ? damping, arH because the representation of the- 
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TABLE 11-1. Values of a and 6 for Three Cases of 
n n 

Unsteady Motions in a Small Motor. 



Cat e (a) 

Case (aa) 

Case (b) 

»1 

- 18. 5 (sec"M 

8. 0 (sec* ^ ) 

-9. 1 


-3o9. 3 

-342.8 

-334. 8 


• 

o 

1 

-583.6 

-566. 5 

“4 

-915. 9 

-889.4 

-871.4 


-12 39.2 

-1262.7 

-1244. 0 

®1 

12. 9 

12. 9 

64. 0 


46. 8 

46. 8 

34. 9 


-29. 3 

-29. 3 

-35. 6 

®4 

-131.0 

-131.0 

-135. 0 


-280. 0 

-280. 0 

-283.0 
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combustion response is the same in the approximate and numerical analy- 
sis, a likely source of the difference in the results is the behavior of the 
nozzle. 

In the numericaJ analysis, the calculations for the entire two-phase 
flow were carried out to the nozzle throat. The influence of the nozzle is 
buried in the results and its contribution to attenuation cannot be determined. 
For the approximate analysis, the influence of the nozzle is represented as 
a surface admittance function evaluated at the nozzle entrance. The re- 
sult (1 1. 1 ) is strictly valid for quasi-steady behavior of a gas only; it has 
been extended to the case of two-phase flow by using the value of y for the 
mixture. The point i 5 that there is reason to expect that the contribution 
of the nozzle is different in the two analyses. 

As a means of comparing the analyses, the value of the attenuation 
constant associated with the nozzle is chosen to obtain he same limiting 
amplitude as found in the numerical analysis. This procedi t:. rests rn 
the observation, elaborated upon further in 5l£., chat the values of the 
limiting amplitudes are quite sensitive to the values of the linear coeffi- 
cients a , 9 
n n 

Not a very large change is required. Equation (11. 1) gives = 

- 153. 25 sec’ If this is changed to - 126. 75 sec" \ then - 8. 00 a id 
the limiting amplitude is about 4. 2 per cent. The values of the 01^ ai d 9^ 
fo; five modes are given in Table 11-1, for the case denoted (aa). Figure 
11-1 shows the amplHudes of the five harmonics < onsidered; the functions 
auid are shown in Figure 11-2 and 11-3, ard a few cycles of the wave 
at limiting amplitude are given in Figure 11-4. 

The attenuation by the nozzle is assumed to vary linearly with the 
ratio J of throat are i to port area, having the value -126- 75 sec"^ hen J 
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has the value (• 132) for case (a). Then with the other required data given 
above» the approximate analysis applied to case (b) produces the results 
shown in the last column of Table 11-1. The initial disturbance is stable* 
the decay constant for the first harmonic being = -9. 12 sec^^. After 
twenty cycles (.02222 sec) the amplitude is roughly 3. 5 per cent according 
to the approximate analysis. The numerical analysis gave an amplitude of 
3. 02 per cent after tv enty cycles. It appeared that the wave may have 
stabilized at a limitirg amplitude, but the calculation was not carried 
further. In view of iiie slow decay found with the approximate analysis, 
it may be that the conclusion based on the numerical analysis, for onlv 
twenty cycles, is incorrect, based on incomplete results. If a true 
limit was indeed reached, then of course the approximate analysis gives 
the wrong result. A limit of 3 per cent would be reached only if a ^ is 
positive, ha\xng a va ue less than 8. 00 . 

That nonlinea." influences are in fact active, even at such relatively 
small aimplitudes, is easily established. Consider a wave having a decay 
constant equal to -9- i2 sec’^ and a frequency of *^00 Hz. According to 
linear behavior, the amplitude would be 4. 08 per cent after twenty cycles 
if the initial amplituc e is 6 per cent. The more rapid decay shown b> the 
nonlinear analysis is evidently due to the transfer* of energy, throuj^^h *:he 
nonlinear processes, from the fundamental oscillation to higher harmonics 
which are then attenuated much more rapidly. Nonlinear particle dariping 
was included in the approximate analysis but, as one should expect fo • the 
small amplitudes arising in these examples, its influence is negligibbr 


small. 
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in 



0.0 0.05 0.10 0.15 0.20 0.25 0.30 


TIME( SEC) 


Figure 11-1, Amplitudes for unstable oscillations in a motor, 
case (aa) of Table 11-1, 





0.0 0.05 0.10 0.15 0.20 0.25 0.30 

’ ,^£(SEC) 


Figure 11-2. The functions A^(t) for case (aa^ of Table 11-1. 




- 0.05 - 0.03 - 0.00 0.02 0.05 
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0.0 0.05 0.10 0.15 0.20 0.25 0.30 


TIME(SEC) 

Figure 11-3. The func tion.s T>j^(t) for c ase (aa) of Table 11-1. 



WfiV 
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11.2 Application to a T-Burner 

The numerical analysis has also been used to analyze a T-burner. 
Both the calculation and the data are for flush cylindrical grains. Data of 
course can be obtained only approximately at that flush condition. 

Most of the properties used are the same as those listed above, with 
the following exceptions. The flame temperature is 3264^K, the particle 
diameter is 3. 0 microns, and the parcimeters in the combustion response 
are A - 8. 8 , B = 0. 67 . 


The T-burner differs from a motor most importantly in two respects: 
the vent is in the centr^r, and there are substantial areas of inert surfcce. 
Apart from possible interactions between the non-.uniform mean flow and 
acoustic velocity, one expects that the center vent should have very little 
influence on the fundamental mode. Indeed, all odd modes, which hav? 
pressure nodes at the center, should not be much affected, indeed, ail 
odd modes, which have pressure nodes at the center, should not be much 
affected. In the calc' lations discussed here, the vent is assumed to have 
no effect on the odd modes. 

Because the even modes have pressure anti -nodes at the center of 
the burner, they will be influenced by the vent, r'or comparison with ':he 
numerical results gi\ en by Levine and Culick (1974), the flow in the vont 
is assumed to be subi^onic and to respond as a plug flow. If the plug has 
effective length L ^ , and the pressure downstream is essentially constant, 
the vent exhausting into a surge tank, then the equation of motion for the 
plug is 


P 


I 

o V 


du’ 

dt 


P' 


(11.7) 


where u' is positive outward from the burner. For harmonic motions, 


wue 
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amplitude of the motion is 


u = - 


p uu-t 

^O V 


( 11 . 8 ) 


,th 


The admittance function for the I mode is therefore given by the formula 

_ f u/li ^ _ i 


A = V ) 

V A ^ / / 


(11.9) 


P''(Po I 

Because the velocity fluctuation of the plug is ninety degrees out of phase 
with the pressure, the admittance function has no real part and there will 
be no contribution to the linear decay constant. The representation of the 
vent as linearized inviscid plug flow introduces no losses. If viscous ef- 
fects or radiation to <he surroundings are included, the phase between the 
velocity and pressur* fluctuations is changed .and there are then associated 
energy losses. 

The influence of the vent is accommodated in the nonlinear formu- 
lation by the term containing R in (3.40), where R is defined by (3. 4<^.). 

Let denote that term so that, if only the effect of the vent is considered, 

th 

the equation for the amplitude "f the n mode is 

.. 2 

T] + (J T1 = F 
n a n v 


^;t it 


P E 

^ o n 


H = -(p +u — j 

V \ o V V — ^ / 

P ^ 


01 . 10 ) 
(U. 11) 

(’1. 12) 


* 


The argument used i: §7 will again be used here >o adapt the admittar :e 
function (11. 9), defined only for harmonic motions, to the general case of 
unsteady motions. Then F^ can be written 

Note the negative sign in R to account for the fact that both u' and u 
are positive outward. 
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F = - ^ fi P Ids . 

V _ 2 dt •) «J n L ■ - ' V n 1(1 v n —2 'n o-l 
E \ n p 3. 

n ’ 

2 

The integrand is nearly constant over the vent, and with L/2 


F = -E 
V nn 


ri -D fi - l'2(u (y-)( ~)a ]n 

n nn n L n^L/^S / v -i n 

c 


(11. 13)' 


th 


According to the rule (6. 10), the linear coefficients for *^he n mode .re 


c 

- X ,S 


^ - -rir' = (f )(s-)\" 


(j 1. 15: 


~n c 

If the mean flow term (proportional to in (11, 14)) is neglected, and the 

admittance function is purely imaginary, as fcr (11.9), then only the • oef- 

ficients 9 for the even modes are affected by the oscillating plug flow, 
n 

It is not diffic'olt to work out the formulas for the linear coefficients 
arising from combusvion, particle/gas interactions, flow -turning, and heat 
losses the lateral boundary. They are given here for the case of c/lin- 
dr’cal grains of length extending from the ends. 

Combustion 

L 






sinZk L, \ 
, , n b^ 

^ ^ -iTTr—J 

n b 


R. 


(r) 


(;1. 16) 


(i) 


9 


n 


R, 


■FT “"n 


(11. 17) 


2 

Here, - ttR^ is tne cross-sectional area of the burner, is the xro3s- 
sectional area of the vent, and is the area of burning surface in one half 
of the burner. 
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Particles 


a 


n 




uu=ou 


n 


(II. 18) 


e 

n 





] 




n 


Flow Turning 

a = zl'Zb)! ^)f, 

n V r / \ L / L 


sin 2k L, 
r b 

2k U 
n b 


(1‘v. 19) 


(J i. 20} 


9=0 

n 


Linear Heat Losses 


( 11 . 21 ) 


a 

n 





x/pt 


yPr 


(II. 22) 


0 = - a (II. 23) 

n n 

Once again, the combustion response is assumed to be given by (7.4) for 
pressure coupling only. 

Both nonlinea particle damping and nonlLiear heat losses are in- 
cluded in the calculations. The first has already been discussed; the formu- 
las are the same as l.aose used in § 8. Nonlinear heat losses are handled 
approximately by using the formulas worked out for a tube, in f*9 , bu' 
weighted by the proportion of lateral surface v/hich is inert. Thus, the co- 
efficient given by (9. -8) is multiplied by (L-2Lj^)/L . Although nonlirear 
heat transfer may be important in setting the vali es of limiting ampli' ades, 
particularly for tests with unmetallized propellants, very little is known at 
the present tune. N< extensive numerical xesults have been obtained spe- 
cifically to determine its influence. 
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Calculations have been done for an example treated numerically by 
Levine and Culick (1974). The propellant is the same as that used in the 
examples discussed in §11. 1. Owing to heat losses, the temperature in the 
chamber is assumed to be somewhat lower, ss3264^K, roughly in accord 
with observati''''.a. Also, the particle diameter is assumed to be 3.0 mi- 
crons, and the parameters in the combustion response are A = 8. 8, B = 

0. 67; these values were chosen somewhat arbitrarily to produce the cor- 
rect growth rate of unstable waves, and should not be taken as representing 
the actual behavior oi the propellant. The computed results will depend 
quite heavily on the combustion response. The main purpose here is to 
compare the approximate and numerical analyses. 

With the data given above, and with L = 23. 5 inches, = 0. 75 inches, 

the values of the constants a and 3 are listed below in Table 11-2 for two 

n n 

values of S. /S . For S, /S = 7.06, two cases are shown: in case (a), 
b c b CO ’ 

there is no influence of the vent; and in case (b), the unsteady flow in the 

vent is treated as a plug flow, as described above. The values of the 

and 0^ are listed in Table 11-2. The values for computed with htat 

transfer neglected agree almost exactly with those deduced from the unstable 

waveforms produced oy the numerical analysis. The amplitudes of the first 

five harmonics, for the case S. /S = 7. 06 and with no influence of the vent, 

D c 

are shown in Figure 11-5. In Figures 11-6 and 11-7 the functions A^ and 
are shown. These should be compared with the results for the example 
for behavior in a motor covered in § 11. 1. For that case, all harmonics 
reached finite limiting amplitudes when the functions A^^ and B^ oscillate 
with fixed amplitudes and frequencies. Here, the first harmonic grows 
without limit, and although the amplitudes of the higher harmonics show a 



RMPLITUDES 

0.12 0.24 0.36 0.48 0.60 
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0.0 0.01 0.02 0.03 0.0« 0.05 0.06 0.07 


TIME(SEC) 


Figure 11-5. Amplitudes for unstable oscillations in a T-burner, 
case (a^ of Table 11-2 (no influence of the vent). 
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TIME(SEC) 


Figure 11-6. The functions A^(t) for case (a) of Table 11-2. 








0.0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 


TIME(SEC) 

Figure 11-7. The fi notions B (t) for case (a) o.< Table 11-2. 
^ n 
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TABLE H-2> 

Values of a and 6 for Two Cases of 
n n 

Unstable Motions in a T- Burner. 



Case (a) 

Case (b) 


73.9 (sec“5 

73.9 


-319.2 

-319.2 


-725.7 

-725.7 

°'4 

-1195.0 

-1195.0 


-1691. 1 

-1691. 1 

0 

”1 

49.2 

49.2 

®2 

7. 5 

-1560.0 

h 

-186.0 

-186.0 

h 

-490. 5 

-1275.0 

h 

-915.0 

-915.0 


decay period after about 0.05 seconds, they too 3ventually grow withort 
limit. It is a curious result that the functions A^(t) and A^Ct) are almost 
identical; note that B^(t) and B 2 (t) are considerably different. Why this 
occurs in this special case is an unanswered question* 

The amplitudes and the functions A^(t), B^^(t) are shown in Figures 
11-8, 11-9, and 11 -lU for the case (b) of Table 11-2. Comparison with the 
results for case (a) shows the influence of the plug flow. As noted above, 
this treatment of the vent provides no losses, but there is a change in the 
relative phases of the harmonics, because the values of 6^ and 0^ are 
quite considerably affected. The most notable feature is the decrease of 
the amplitude of the first harmonic in the period . 05 - . 09 sec. , after which 
it again grows without limit. The values are unrealistically h?gh, much 




AMPLITUDES 

0.14 0.28 0 



0.0 0.02 0.04 0.06 0.08 O.IO 0.12 

TIME( SEC) 


Figure 11-8. Amplitudes for unstable oscillations in a T -burner, 
case (b) of Table 11-2 (unsteady fTw in the vent 
treated as a plug flow). 
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0.0 0.02 0.04 0.06 0.08 0.10 0.12 

TIME(SEC) 


Figure 11-9. The functions A^(t) for case (b) of Table 11-2. 
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TIME(SEC) 


Figure 11-10, The Functions B^(t) for Case (b) of Table 11-2. 





- 99 . 


larger than those obs^rveo for the tests which were the basis of these cal- 
culations* 

It is a puzzling difficulty at the present time that for apparentl> 
realistic values of the data required^ the calculated amplitudes for motions 
in a T -burner are much too large and may even grow without limiting. This 
has been the case with other examples not included here. A similar diffi^ 
culty was encounterer: in the numerical calculations. For the case (b), 
the numerical results seemed to indicate that the amplitude was limiting, 
at a value of roughly . 37, at approximately . 08 seconds; after this time, 
the amplitude of the fundamental mode according to Figure 11-8 begins to 
grow again. The numerical results were not carried out further. 

Thus, for the examples treated here, limiting behavior has not been 
obtained for the unstable motions in a T -burner. Even with the growth con- 
stant for the fundamental mode arbitrarily reduced, the amplitude grows 
without limit; as ttj, tends to zero, a longer time is required for the ampli- 
tude to reach unity. This is a striking contrast with the results obtained 
for a motor. In both cases, reasonable vadues of the required data have been 
used. But as shown by Tables 11-1 and 11-2, the relative values of the a^, 

0^ for the various harmonics are considerably different. This is clearly 
the origin of the difference in nonlinear behavior, but the details are not 
known, and no generalizations can be offered. 

The difference between cases of limited and unlimited growth appears 

strikingly in the functions A (t) and B (t). K the value of 0 is non-zero, 

^ ^ n n n 

both functions oscillate (see §12.2) when the amplitudes limit. Otherwise, 
the behavior has no covious pattern. The contrast may be seen by comparing 
the results of §11.1 with those shown here. This question is remarked upon 
further in 512.2, but a general auiswer is yet to be found. 
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XII. THE CONNECTION WITH LINEAR STABILITY ANALYSIS AND 
THE BEHAVIOR OF APPROXIMATE NONLINEAR SOLUTIONS 

Most problems of combustion instability in solid propellant rockets, 
liquid propellant rockets, and in thrust augmentors have been studied by 
using linear stability analysis. It has been emphasized in the formal con- 
struction of the analysis given here* and in the examples covered in §s 8 
and 11, that the linear coefficient is exactly the familiar decay or 
growth constant for the n^^ harmonic. One of the appealing features of the 
results (cf. eqns. (4. 18), (4. 19), and (6. 10) -(6. 12)) is that the corre- 
spondence is so readiJy established and obvious. In practice, this means 
that the numerical results obtain »d from linear stability analysis may be 
used directly in the approximate nonlinear analysis. Indeed, for the ex- 
amples given in § 11, the cost of doing the necessary linear calculations 
was comparable to that of the numerical solution to the nonlinear equations. 
In § 12. 1, the relation between the present analysis and the more familiar 
linear stability analysis is elaborated upon further. 

One of the interesting questions which arises in connection with the 
nonlinear analysis is; under what conditions will the motions reach a 
limiting simplitude? Put ainother way, for what rcdinges of the linear coeffi- 
cients {a^ , 0^3 will the system of coupled nonlinear oscillators execute a 
limit cycle? / ^ the 'iresent time, the answer to this question is not known, 
and therefore it is also not possible to offer any i<eneralizations about the 
values of limiting amplitude, and how they depend on the linear coefficients. 
A few remarks on this subject, which merits fur her work, are given in 
5 12 . 2 . 
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12. 1 The Connection with Linear Stability Analysis 

All analyses of linear stability are ultimately related in one way or 
other to the inhomogeneous wave equation with inhomogeneous boundary 
conditions, having the form 


a 8t^ 

-h 

(12. 1) 

n . Vp' = -f 

• 

(12.2) 

For the formulation given earlier. 

h stands for the sum of h 

U 

and the 


linear part of h^ in eq. (3. 8); f represents f^ , f^ and the linear part of 
f in (3. 9,.. 

V 


The linear stc.bility of normal modes is studied by assuming tliat 
p' has the form 


p’(?.t) = ' 


(12.3) 


and the problem comes down to determining the mode shape p(r) and the 
complex wave number k : 

k = ~ (ou-ia^) (12.4) 


a 

where is the grov'th constant for the perturbed n^^ mode. Because h 
and f are linear in perturbations (they both contain the velocity u' well 
as p'), they may be written in the form 


h(r, t) = n(r)e , f(r, t) = f e 

(12. 5) 

Then (12, 1) and (12,2) become 


V p -1- k p = h 

(12. 6a) 

:l • V p = -f 

(12. 6b) 


A A 

The functions h and f are of first order in the Mach number of the mean 

flow, assumed to be small. It is then a relatively simple matter to deduce 
2 

a formula for k , valid to first order in the Mach number (Culick 1975, 
for example): 
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c 

n 

Equation (12. 7) is based on an iterative procedure; this first approximation 
requires substitution of the unperturbed mode shape acoustic 

^ A 

quantities appearing in h and f . 

Because ^ 1 « the real eind imaginary parts of (12. 7) can oe 


written 


2 

RefJf^hdV + ^ t|;^f dS} 


n 


p E 
^o n 


a 


n 


2p ujE 
^o n 


2 - Jm{Ji^^hdV + |Jl|(^fdS} . 


( 12 . 8 ) 


(12. V) 


It is helpful to make ’-he discussion more specific by considering only one 
term for a specific problem. For this purpose, the contribution from 
surface combustion i-i a cylindrical motor serves quite nicely. Then ac- 
cording to eq. (3. 4), 

t = pQ 15t~ ’ ^ ’ (12.10) 


so 

A 

f = iakp^u* n . (12. 11) 

Then with the definitions introduced in § 7.6, and if non-isentropic tem- 
perature fluctuations are ignored, (12, 11) can be expressed in terms of 
the response function: 

t" = . (12.12) 

Substitution into (12. 8) and (12. 9) gives 

= yUj^(^)(2(i)^)Rj^^^^ , 

c 


(12. 13) 
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- 'YUj^( )R^ (12. 14) 

C 

Now because it is a requirement, in order for th^ • calculation to be valid, 
that the perturbations be small, ou must not differ much from , so 
(12. 13) can be written 

6it)^ = ^ ^ - (12. 15) 

c 


More generally, (12. 8) gives 

2 2 

cu -(W 


-e = 


n 


ToT 


w (JU-(J0 


~2 

a 


n 


n 


2p uu E 
^o n n 


^Re{Ji|;^hdV + |^i|(^fdS} . (12.16) 


For most work in the past dealing with combustion instability, the 
frequency shift due to perturbations has been ignored. It is normally too 
small to be distinguished in experimental work. However, it happens, as 
the example in § 11.2 suggests, that the nonlinear behavior -- in particular 
the limiting amplitude -- be quite sensitive to those characteristics 
which cause a linear frequency shift. The reason is that the phase rela- 
tionships between the various "armonics are aiffected. 

The formal connection between the linear stability analysis and the 
nonlinear analysis based on the method of averaging may be seen most 
quickly by comparing eqs. (12. l4) and (12. 13) or (12. 15) with eqs. (7. 9) 
and (7. 10), which for a cylindrical grain become (11.4) and (11. 5). For 


this special case, the values of are the same, and 


6^ = . (12.17) 

n n 

The fact that the linear coefficient 0^ is the negative of th^ frequency shift 

can be shown quite easily by starting with the definition of ri^ : 

r\ = G 8in(uu 4) ) = A sin{(j) t) + B cos('ju t) . (12. 18) 

nnnn n n n n 

The frequency of the actual wave (i. e. , when perturbed from the mode 



-104. 


RETBODUCroilJIV 

PAe* B POOR 


having frequency ir^ ) is 

d 

.; = •^r (uu t + «j> ) = U) + i , 
dt ' n ^n' n ^n ’ 

so 


6(1) = i 

n ^n 


(12. 19) 


Now tan 4> = B /A , and differentiation with respect to time gives 
n n n o 


B / B \ A 

1 _ n f n n\ n 

*^n "" A_ B " A . 2^ 2 

n n n A 


( 12 . 20 ) 


The linear parts of ths equations for A^ and B^ are, from (4. 18) and 


(4. 19), 


dA 

: a A + 9 B 
dt n n n n 


dB_ 


n 




= a B - 9 A 
n n n n 


(12.; la) 
(12.21b) 


From these it follows that 


A ^+B ^ 


n '“n _ ft f n n ') 

■ A * "®n \ A„B / ’ 


n "n n n 

and substitution into '.12. 20) gives the result 


6u; = = -e^ . (12.22) 

21 n n 

This relation has been referred to in § 8; numerical results for the shift 
of frequency (or dispersion) associated with gas/particle interactions are 
given in Table 8-1. 

The main point here is that numerical results obtained with the ap- 
proximate nonlinear analysis have shown that the nonlinear behavior i.iay be 
quite sensitive to the quantity 9^ = f which has largely been ignored 

in linear stability analysis. 

It should be clear that the relation 6 = -6(ti is true in general, 

n n 

and not just a specia* result for surface combustion. This can be seen di- 
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rectly by returning to eq. (3. 15) for , and retaining only the linear 


terms incorporated in h and f used above: 


v%\ = 


-2 

a 


n 


rt I ♦„Mv + . 


(12.23) 


For harmonic motions, ^ = 'Hjj exp(iakt) and because both h and f are 
proportional to , the factor can be dropped, so that (12. 18) becomes 


-(ID -ia)^+(D ^ 
n 




^o n 


2-{Ji^hdV+ ffd^fdS] . 


The real and imaginary parts are once again (12. 8) and (12. 9) 

When the method of averaging is applied to (12. 18), the right hand 
side is expressed as a combination of terms depending on r\. or T). , eqns. 
(4. 1) and (4. 2); with only the linear terms retained, (12. 18) is 


- £ OD .T). + E .T1.] . 


(12.24) 


n+(DTi=-D‘n-E‘n _ 

n n n nn n nn n ni i ni 1' 

The analysis developed in §4 leads to the first order differential equations 

(4. 18) and (4. 19) for and . Those equations show that in general 

there are terms representing linear coupling betv\-een the modes providing. 


of course, that the coefficients E ., E . don't vanish. These terms do not 

m ni 

arise in the usual linear stability analysis because the restriction is en- 
forced from the beginning that the motion consists only of a wave having a 
single frequency. If all the and E^ are zero, or for the case of 

longitudinal modes, the equations for the and are eventually found 


to be (12.21a, b), with 

a 


nn 


nil 


n 


(12.25) 


The only remaining point to establish is that the and 6^ are the 
same as defined by eqns. (12. 9) and (12. 16). This is easily shown in 
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special cases. To see that the equalities are true in general, begin with 
the two expressions <12. 23) and (12. 24) for the right hand sides, and 


D . = E . = 0 : 
ni ni 


D +E Y! ^ 
nn n nn n 


-2 

a 


p E 
*^o n 


j{r»__hdvt jr,__fds} . 


(12.26) 


Now the same argument is used for the left hand side of (12.26) as that 
introduced first in the paragraph preceding eq. (7. 7). Namely, the ap- 
proximation is made tor harmonic motions. ^ ^^n^n * another way, 

because the D , E are already small quantities, and real, then for 
nn nx 


harmonic motions. 


D r? 
nn n 


vx D r\ e 
n nn n 


ix t 
n 


The real and imaginary parts of (12. 26) are then 


D 

nn 


~2 

a 


p X' £ 

^ ^ n n 


2-<inCj<r^hdV+ IJu/dS] 


-2 


E « C '*v hdV + rr t fd£} 

nn _2 '-•'’n JJ n-* 

P E 


(12.27) 


(12. .18) 


n 


Comparison with (12. 8) and (12. 16) shows immediately that 


and 6 


-E /2x a^ required, 
nn n ^ 


12. 2 Remarks on thr Behavior of Solutions Obtained with the Approximate 
Nonlinear Analvsis 


For the example given in §11. in which aJl amplitudes reach 
limiting values, the functions A^(f) snd both oscillate, with the same 

amplitude and a phase difference of tt/2 after the limit condition has been 
reached. This sort >f behavior has been most commonly found in cases 
for which the motion eventually executes a limit cycle. There are examples, 


as shown below, in vhich the A and B as well as the amplitudes 

n n ^ 
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vA +B tend in the imiit to constant values. If a limit cycle does not 
exist, there seems to be a wide variety of possible behavior. One example 
has been given in §11.2; several more will be included in this section. 

According to Lie discussion in §12. 1^ the linear behavior of the 


functions A^(t) and t>^(t) is oscillatory with exponential growth or dev^ay. 
The frequency of the oscillation is 0^, the negative shift of the normal 
frequency To see this explicitly, note that b/ definition. 


A (t) = G cos 4> 
n n n 


B (t) - G sin 4> 
n n ^n 


(12. ?9) 


According to (12. 22), a = -0 , so <{> is a constant minus 0 t. Hence, 

n n n n 

A^ and B^ both oscillate at the frequency 0^. Another way of establishing 
this result, and at th#-- same time incorporating the nonlinear inruences, 
may be had by combining the equations for A^ and B^; write these in the 
form 


= a a + 0B +f 

dt n n n n n 


•’.0) 


= CL B - 0B +6 
dt n n n n ®n 


( 12 . 51 ) 


where stand for nonlinear terms and, in the general case, for linear 

coupling terms as well. 

Differentiate (12. 30) with respect to time and then use the two equa- 
tions to eliminate B^ and B^ . A similar procedure Ceui be applied to 
(12.31). The two second-order equations for and B^ are produced: 


— -2a + (0 ^+0 ^)A = f -a f +0 g , 

n dt n n ' n n n n n®r. 


( 12 . ^ 2 ) 


2 

d B d B 7 7 

_!1 -2a^-^ +(a^+0/)B^ = L-O'ngn-Vn 

ndt nnn nnnnn 
dt 


( 12 . ^ 3 ) 
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These equations represent two coupled nonlinear oscillators asso- 
ciated with each of the normal modes. The homogeneous solutions for the 
lins ir operators are 

a t iS t 

- e “ e " , (12. 34) 

e.o hat for an unstable mode, the A^, exhibit oscillatory growth or de- 
cay. For cases in which the amplitudes A^, B^ show the limiting behavior 
note 1 above, the nonlinear terms on the right hand side must be taken Into 
acct'unt. If they have relatively small influence oa the frequency, then in 
the limit cycle, both and B^ oscillate with frequency 0^ , the neg itive 
of tht; frequency shift of the frequency for the primary oscilla.ors repre- 
sented by the r, . 

' 'u 

But if the motion does go into a limit cycle, the nonlinear ternr 3 on 
the right hand sides ci (12. 32) and (12. 33) must have some influence. A 
pu 'ely oscillatory behavior would be produced by (12. 32) if the term 

is exact^> compensated and if the remainder of the right hand side 
cortrioutes a term proportional to A^. It would be interesting and quite 
likely very useful to determine the general conditions under which solutions 
of that type do or do not exist. 

Etiuation ^12. ^4} shows that if 6 =: 0 . the functions A and B do 

n n n 

not oscillate w 'r'n the motion is linear. There are cases in which tha^. be- 
havior persists when nonlinear influences are ac< ounted for (see case (ii) 

below); o.nd there are conditions under which the A and B are not rarely 

n n 

os- iatory even though 0^ ^ 0 (case (iii) below). Again, quantitative rjener- 
ctlizations would probably be very useful for practical applications. 

Cne way of attacking this problem would to use the method nf 
avera^;,ing to solve (1 .32) and (12. 33). This would not give general r esults 
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covering all possibilities, but it might yield information at least in respect 

to the conditions for oscillatory behavior of A^(t) and B^(t) . 

No formal results have been obtained for the general behavior of 

the solutions. The practical importance of pursuing the problem lies in 

the possibility for defining the ranges of values of the linear coefficients, 

0.^ and 6^ , for which limit cycles exist; and for determining the dependence 

of the limit amplitudes on those parameters. 

In the absence of general results, the problem has been studied in 

a very modest way simply by carrying out calculations for special cases. 

Seven examples are included here; the values of a and 6 are listed in 

n n 

Table 12-1. Case (i) has been chosen as the reference because it exhibits 

the relatively nice oscillatory behavior of the functions A^(t) and B^(t) . 

For each case, die results are shown for the anplitudes and for the 

the B (t) behave qualitatively like the A (t) in all cases. The calculations 
n n 

were done for five modes, but to make the figures a bit less crowded, only 
the curves for the first three modes are shown ic Figures 12-4 to 12-7. 

Note that in tlie reference case, Figure 12-1, the first mode io un- 
stable and all higher modes are rather heavily damped with = -IOC sec“^. 
Case (ii). Figure 12-2, shows the change to non-oscillatory growth ol the 
An(t) when = 0 • Cases (iii) through (vi). Figures 12-3 to 12-6, show 
the influence of makmg each of the higher modes neutrally stable. The re- 
sults reflect at least partly the strength of the coupling between the £i-*st 
and each of the higher modes. Evidently the higler odd modes are m >re 
strongly coupled to the fundamental than are the higher even modes* \Vith 
a i or equal to zero, the motions grow without limit; but when or 
are zero, the mo-ion does reach limit cycles, although they differ in 
detail. Even the qualitative aspects are not obvious from the equations 
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(6. 16) - (6.20) which have been solved. 

Figure 12-7 for case (vii) is included as an illustration of somewhat 
more erratic behavior which often ensues if more than one mode is un- 
stable. Still more extreme cases have been encountered, but nothing can 
be accomplished by mustering all the results which have been obtainec. 



TABLE 121-1. Values of a and e for j''igures 12- 
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AMPLITUDES 

0.01 0.02 0.03 0.04 0.05 
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Fig. 
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12 -2(b) The fuiif. tions A^(t) for case (ii)# 
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Fig. l2-4(a) AmpUtides for case (iv); values of 6^ and (X..^ 
same ts case (i), except = 0 . 
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Fig, l2-4(b) The fan :tions A^(t) for case (iv). 
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Fig. 12.5(a) Amplitudes for case (v); vadues of 9 

n 

as for case (i) except = 0 . 


and a 

n 


same 
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Fig. 12.5(b) The furctions A^(t) for case (vi). 
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Fig. 12 -6(a) Amplitudes for case (vi); values of 6^ and same 
as case (i) except = 0 . 
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Fig* l2-6(b) The functions for case (vi). 
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Fig. l2-7{a) Amplitudes for case (vii); values of 9^ and same 
as for case (i) except = 0 and = 16 . 
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Fig. 12-7(b) The Functions for case (vii). 
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